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It has been known since the pioneering work of Bernal, Fowler and Pauling that common, hexagonal (Ih) 
water ice is the archetype of a frustrated material : a proton-bonded network in which protons satisfy strong 
local constraints — the “ice rules” — but do not order. While this proton disorder is well established, there is 
now a growing body of evidence that quantum effects may also have a role to play in the physics of ice at low 
temperatures. In this Article we use a combination of numerical and analytic techniques to explore the nature of 
proton correlations in both classical and quantum models of ice Ih. In the case of classical ice Ih, we find that 
the ice rules have two, distinct, consequences for scattering experiments — singular “pinch points”, reflecting 
a zero-divergence condition on the uniform polarization of the crystal, and broad, asymmetric features, coming 
from its staggered polarisation. In the case of the quantum model, we find that the collective quantum tunnelling 
of groups of protons can convert states obeying the ice rules into a quantum liquid, whose excitations are 
birefringent, emergent photons. We make explicit predictions for scattering experiments on both classical and 
quantum ice Ih, and show how the quantum theory can explain the “wings” of incoherent inelastic scattering 
observed in recent neutron scattering experiments [Bove et al, Phys. Rev. Lett. 103, 165901 (2009)]. These 
results raise the intriguing possibility that the protons in ice Ih could form a quantum liquid at low temperatures, 
in which protons are not merely disordered, but continually fiuctuate between different configurations obeying 
the ice rules. 


I. INTRODUCTION 

We learn as children that matter can exist in three differ¬ 
ent phases — solid, liquid, and gas. This concept is usually 
introduced through the example of water, familiar as a liquid 
(water), a gas (steam) and a solid (ice). However, at least as 
far as its solid phase is concerned, water is a spectacularly 
unusual material. At atmospheric pressure, water molecules 
freeze into a structure known as “ice Ih”, illustrated in Fig.[^ 
Ice Ih is remarkable in that the oxygen ions (0^“) form an 
ordered lattice, while the protons (H+) lack any kind of long- 
range order — in flat contradiction with the usual paradigm 
for solids. 

This extraordinary property of water ice was first elucidated 
more than 80 years ago, by Bernal and Fowler m Bemal 
and Fowler argued that ice should be viewed as a molecular 
solid, in which distinct water molecules are bound together 
by hydrogen bonds. Each water molecule forms four such 
hydrogen bonds, and as a result, the proton configurations 
obey strong local constraints, commonly referred to as the 'fee 
rules'' mm. The ice rules lead to strong correlations between 
protons, but can be satisified by an exponentially large num¬ 
ber of different proton configurations mm. As a result, the 
protons remain disordered, and possess an extensive residual 
entropy. This “ice entropy” is observed in experiments on Ih 


water ice and persists down to the lowest temperatures 
measured, in apparent defiance of the laws of thermodynam¬ 
ics. Eighty years on, these striking discoveries continue to 
exert a profound influence on research into water ice O, and 
a wide range of other materials IbHTTl. 

Recent experiments by Bove et al CD suggest a new twist 
on the behaviour of protons in ice Ih — not only are protons 
disordered, but they remain mobile, even at temperatures as 
low as 5 K. This might seem surprising, since any attempt to 
move a proton will lead to a violation of the ice rules, at con¬ 
siderable cost in energy El. This problem is avoided, how¬ 
ever, if the proton dynamics consists of coherent collective 
quantum tunnelling on hexagonal plaquettes, of the type il¬ 
lustrated in Eig. 1^ This mechanism for proton dynamics in 
ice Ih flnds support from ah initio calculations |[T9|[20l, with 
recent results suggesting that, while the high-temperature dy¬ 
namics proceeds via single-proton hopping, collective motion 
around loops becomes important at low temperatures 1201 . 
And an analogous correlated tunnelling of protons has re¬ 
cently been observed in an artifleial assembly of four water 
molecules ca. 

There are, in fact, many classes of system whose low- 
temperature physics is subject to strong local constraints, sim¬ 
ilar to those found in water ice. The most celebrated of these 
are the magnetic systems known as the "spin ices", whose 
low-temperature spin conflgurations are in correspondence 
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FIG. 1. (Color online). Crystal structure of hexagonal (Ih) water ice. Water ice can be viewed as a frozen configuration of water molecules, 
satisfying the Bemal-Fowler '‘ice rules” (BID. in which each oxygen (red sphere) forms two short, covalent bonds, and two long, hydrogen 
bonds with neighbouring protons (white spheres). Oxygen atoms form an ordered lattice, belonging to the hexagonal space group PGa/mmc, 
with a 4-site primitive unit cell. Protons do not show any long-range order, (a) Structure viewed perpendicular to the hexagonal symmetry 
axis (the crystallographic c-axis). (b) Structure viewed along the hexagonal symmetry axis. 
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(a) Tunnelling on a plaquette within hexagonal plane. 





(b) Tunnelling on a plaquette connecting hexagonal planes. 

FIG. 2. (Color online). Collective quantum tunnelling between different proton configurations satisfying the ice rules in hexagonal (Ih) 
water ice. The crystal structure of ice Ih contains two distinct types of hexagonal plaquette, each containing six oxygen atoms (red spheres). 
Tunnelling between different ice configurations is mediated the by correlated hopping of protons (white spheres) around such a plaquette. 
Where protons form an alternating sequence of short (S) and long (L) bonds (e.g. S-L-S-L... -S-L), it is possible to tunnel to another, 
degenerate, ice configuration in which long and short bonds are interchanged, (i.e. L-S-L-S.. .-L-S) |[T8H20l . (a) Tunnnelling on hexagonal 
plaquette of “type I”, in th^lane perpendicular to the hexagonal symmetry axis. This process has a matrix element gi in our model for proton 
dynamics, Tunnnelling on hexagonal plaquette of “type IF, connecting different hexagonal planes. This process has a 

matrix element g 2 . 
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with the proton configurations in water ice Ii|7l[l0|. Ice¬ 
like physics also arises in models of frustrated charge order 
Gil US, proton-bonded (anti-)ferroelectrics |[T^fT5l[2TH25l . 
dense polymer melts 1^ and dimer models (271 [28l. In these 
systems, violations of the ice rules take on the character of 
fractionalised charges This point has attracted particu¬ 
lar attention in the case of the spin ices, since these charges 
behave as effective magnetic monopoles (29^3^ . And, in 
fact, a direct analogy may be drawn between these magnetic 
monopoles and ionic defects in water ice (^13311^ . 

The effect of quantum tunnelling, of the type proposed by 
Bove et al (m, has been studied in a range of other ice¬ 
like problems, with striking results. Quantum tunnelling has 
been shown to give rise to a quantum “spin-liquid”, compris¬ 
ing a coherent superposition of an exponentially large num¬ 
ber of states obeying the ice rules, in models derived from 
spin ice EUBD. Equivalent quantum liquids have also been 
found in three-dimensional quantum dimer models Il42l - l44ll . 
In both cases, the low-energy excitations of the liquid are 
gapped, fractionalised, “charges” and gapless, linearly dis¬ 
persing, “photons”, in direct correspondence with the theory 
of electromagnetism |[^l42ll . 

A similar “electromagnetic” scenario has also been dis¬ 
cussed in the context of a simplified model of water ice by 
Castro Neto et al. (451 . Finding that quantum fiuctuations 
drive the protons to order in a two-dimensional model of 
water ice C51I76M1, these authors argued, by extension, 
that quantum effects in a three-dimensional water ice could 
drive a finite-temperature phase transition between a low- 
temperature proton-ordered phase, and a high-temperature 
proton disordered phase. Within this scenario, ordered and 
disordered phases of protons in hexagonal water ice — ice Ih 
and ice XI a — would correspond to the confined, and de- 
confined phases of a compact, frustrated U{1) lattice-gauge 
theory (46l|47l, similar to that proposed in the context of quan¬ 
tum spin ice isnsi. 

In the work which we present here we do not attempt to 
establish the conditions under which the protons in hexago¬ 
nal water ice order, but instead seek to characterise their de- 
confined, disordered phase. To this end we develop theories 
which describe the proton correlations in both classical and 
quantum models of hexagonal (Ih) water ice, making explicit 
predictions for scattering experiments. 

In the case of classical ice Ih, we find that the ice rules have 
two, distinct, consequences for proton correlations, directly 
visible in scattering experiments. Firstly, algebraic correla¬ 
tions of the uniform polarisation lead to “pinch points” — sin¬ 
gular features in scattering — visible in a subset of Brillouin 
zones. Secondly, exponential correlations of the staggered po¬ 
larisation lead to broad, assymetric features in a different sub¬ 
set of Brillouin zones. This analysis provides new insight into 
diffuse scattering experiments on ice Ih (48l[49l, and makes 
explicit the differences between ice Ih and spin ice, or (cubic) 
ice Ic. It also provides the starting point needed to construct a 
theory of ice Ih in the presence of quantum tunnelling. 

In the case of quantum ice Ih, we find that it is possible 
to describe the proton configurations in terms of a lattice- 
gauge theory in which the connection with electromagnetism. 


long-implied by the ice rules, is made explicit. The decon- 
fined, proton-liquid phase of this theory is shown to support 
two types of excitation — gapless, linearly-dispersing pho¬ 
tons with birefringent character, and weakly-dispersing opti¬ 
cal modes, corresponding to local fiuctuations of the electric 
polarisation. The predictions of this lattice gauge theory are 
shown to be in quantitative agreement with the results of vari¬ 
ational quantum Monte Carlo simulations. 

Throughout our analysis we will emphasise the ways in 
which proton-correlations in hexagonal water ice differ from 
spin-correlations in cubic spin ice, paying particular attention 
to the new structure which arises in both the classical contin¬ 
uum field theory, and the quantum lattice-gauge theory. 

The remainder of the Article is structured as follows : 

In Section ini we introduce a model of ice Ih which includes 
tunnelling between different proton configurations obeying 
the “ice rules”, and describe those aspects of the symmetry 
and geometry of the ice Ih lattice which are important for our 
discussion. 

In Section we develop a coarse-grained, classical field 
theory describing the correlations of protons in ice Ih, in the 
absence of quantum tunnelling, and compare this with the re¬ 
sults of an equivalent, lattice-based calculation, developed in 
Appendix The two approaches are shown to agree in the 
long-wavelength limit. The lattice-based calculation is also 
used to make explicit predictions for scattering experiments 
on a ice Ih. These are summarised in Fig.[^ 

In Section IV we develop a quantum (7(1) lattice gauge 
theory describing the correlations of protons in ice Ih, in the 
presence of quantum tunnelling. We construct the excita¬ 
tions of this lattice gauge theory, which include birefringent 
emergent photons, and make quantitative predictions for their 
signature in inelastic scattering experiments. Key results are 
summarised in Fig.[7]and Fig.[^ 

In Section [V| we discuss our results in the context of pub¬ 
lished experiments on Ih water ice. We consider in particular 
the long-wavelength features seen in diffuse neutron scatter¬ 
ing, the incoherent inelastic neutron-scattering experiments 
of Bove et al ca, and the thermodynamic properties of ice 
at low temperatures. 

We conclude in Section Ivg with a summary of our results 
and a brief discussion of open questions. 

Technical aspects of this work, including comparison with 
classical and quantum Monte Carlo simulations, are devel¬ 
oped in a series of Appendices. 

Appendix [A| provides details of the ice Ih lattice and of the 
coordinate system used to describe it. 

Appendix [^provides details of the derivation of the classi¬ 
cal, continuum field-theory introduced in Section [nl[ 

Appendix provides details of the microscopic calcula¬ 
tion of the proton-proton correlations in a classical ice Ih, 
described in Section 11111 

Appendix provides explicit relationships between the 
measurable proton correlations and the correlations of the 
fields appearing in our continuum theory. 

Appendix [^provides details of the quantum (7(1) lattice- 


gauge theory introduced in Section IV 


Appendix provides a comparison between the analytic 
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FIG. 3. (Color online). Different types of bond within the ordered oxygen lattice of hexagonal (Ih) water ice. The Ih lattice can be viewed 
as a stack of buckled, honeycomb lattices, composed of bonds which are symmetric under inversion about the centre of the bond. These 
honeycomb layers are connected by bonds running parallel to the hexagonal symmetry axis, which are symmetric under reflection in the 
plane perpendicular to the bond. This lattice is bipartite — i.e. it may be divided into two sublattices — here colored red and blue, (a) 
Centre-symmetric bond, viewed from a direction perpendicular to the bond, (b) Centre-symmetric bond, viewed along the bond direction, (c) 
Mirror-symmetric bond, viewed from a direction perpendicular to the bond, (d) Mirror-symmetric bond, viewed along the bond direction. 


predictions for proton correlations described in Section[nI|and 
Section|IV| with Monte Carlo simulation of both classical and 
quantum ice Ih. 

Appendix provides the derivation of the structure fac¬ 
tor for inelastic, incoherent neutron scattering as measured by 
Bove et al in ca. 


II. CLASSICAL AND QUANTUM MODELS OF 
HEXAGONAL (IH) WATER ICE 

The key to understanding the structure of water ice is the 
realization, due to Bernal and Fowler, that water molecules 
retain their integrity on freezing (Q. It follows that each oxy¬ 
gen remains covalently bonded to two protons, while at the 
same time forming two, weaker, hydrogen bonds with protons 
on neighbouring water molecules. In the frozen state, the oxy¬ 
gen atoms form an ordered lattice, held together by interme¬ 
diate protons, each of which forms one long (hydrogen) and 
one short (covalent) bond with a neighbouring oxygen. This 
type of bonding favours a tetrahedral coordination of oxygen 
atoms, but does not select any one structure, with 17 different 
forms of ice crystal known to exist la. 

Hexagonal (Ih) water ice, illustrated in Fig. is the most 
common form of water ice, formed at ambient pressure. In it, 
oxygen atoms form a crystal with the hexagonal space group 
P 63 / mmc. This structure can be thought of as a set of buck¬ 
led honeycomb lattices, composed of centre-symmetric bonds 
of the type shown in Fig. |^a) and (b). These honeycomb 
layers are linked by mirror-symmetric bonds, parallel to the 
hexagonal symmetry axis, of the the type shown in Fig. |^c) 
and (d). The two types of bond have almost exactly the same 
length, leading to a near-perfect tetrahedral coordination of 
oxygen atoms ilSl ISOUS^ . The primitive unit cell of ice Ih con¬ 
tains 4 oxygen atoms or, equivalently, 4 water molecules, with 
8 associated protons. This should be contrasted with the 2- 
site primitive unit cell needed to describe the diamond lattice 
of oxygen atoms in cubic (Ic) water ice (space group Fd3m), 


or its magnetic analogue, spin-ice El. 

While the oxygen atoms in ice Ih form an ordered lattice, 
protons do not. The way in which water molecules bond to¬ 
gether does not select any one proton configuration m, but 
rather an exponentially large set of 



proton configurations, where A" is the number of oxygen 
atoms (equivalently, water molecules) in the lattice OH. As 
a consequence, the protons do not show any long-range order. 
While unusual, extensive degeneracies of this type are by no 
means unique to water ice, occurring in “spin ice” Clllol, 
problems of frustrated charge order on the pyrochlore lat¬ 
tice lEiini, and in a wide range of problems involving the 
hard-core dimer-coverings of two- or three-dimensional lat¬ 
tices EiiEa. 

The principles governing the arrangement of protons in wa¬ 
ter ice are neatly summarised in the Bernal-Fowler 'Hce rules'' 

EElEl : 

1. Each bond between oxygen atoms contains exactly one 
proton. 

2. Each oxygen has exactly two protons adjacent to it. 

Since the ice rules define one of the simplest models with an 
extensive ground-state entropy, they have proved a rich source 
of inspiration for statistical studies EEIllsa, particularly in 
two dimensions, where the corresponding “six-vertex model” 
can be solved exactly 156^621 . 

Violations of the ice rules cost finite energy, and fall into 
two types. Violations of the first ice rule, double-loaded or 
empty bonds, are known as Bjerrum defects |[63]| . Violations 
of the second rule occur where a proton is transferred from one 
water moleule to another, creating a pair of “ionic” defects, 
hydroxil (OH“) and hydronium (H 3 O+) d O |33l. These 
ionic defects are in direct correspondence with the fractional 
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charges found in models of frustrated charge order on the py- 
rochlore lattice ca, and the magnetic monopoles observed in 
spin ice 129^321 . 

While the proton configurations which satisfy the ice rules 
do not exhibit long-range order, they do possess a definite 
topological structure. This, and the connection between the 
ice rules and electromagnetism, is most easily understood if 
proton configurations are viewed in terms of a conserved fiux. 

The mapping onto a flux representation starts with the ob¬ 
servation that, in the absence of Bjerrum defects, all proton 
configurations can be represented by a set of arrows on the 
bonds of the oxygen lattice. Each arrow points in the direc¬ 
tion of the displacement of the proton from the midpoint of 
that bond, as is illustrated in Fig. Since the oxygen lattice 
is bipartite (i.e. it may be divided into two sublattices), each 
arrow can be thought of as the fiux of a vector field P from one 
sublattice to the other. The second ice rule then amounts to the 
condition that fiux is conserved, i.e. there are two incoming 
arrows and two outgoing arrows at each oxygen vertex. This 
in turn can be viewed as a zero-divergence condition on the 
fiux, 

V • P = 0 , ( 2 ) 

true at every vertex of the lattice. 

The fiux representation of water ice and related systems is 
an approach with a long history IT^ . and is particularly useful 
in two dimensions, where ice states map onto the exactly- 
soluble six-vertex model || 6 ^ . The existence of a zero- 
divergence condition, Eq. Q, suggests a natural analogy with 
electromagnetism, which we will explore further in the re¬ 
mainder of this article. The fiux representation also plays 
a crucial role in understanding scattering experiments in, 
since both the electric polarization of a given bond, and the 
distribution of mass on that bond, are determined uniquely by 
the fiux of P. 

The topological structure of the ice states becomes evident 
once periodic boundary conditions are imposed on the lattice. 
In this case, the local conservation of fiux (second ice rule) 
gives rise to a distinct set of global topological sectors, with 
definite, quantised, flux through the (periodic) boundaries of 
the crystal. In a real crystal, with open boundary conditions, 
these topological sectors correspond to the different possible 
values of the electrical polarisation of the crystal. The very 
high dielectric constant of water ice 0 , can therefore be in¬ 
terpreted as evidence of fluctuations between different topo¬ 
logical sectors ll^[3^l64|[65]| . 

Quantum mechanics enters into the physics of ice through 
the (quantum) tunnelling of protons from one configuration 
to another. Since it is energetically expensive to violate the 
ice rules, tunnelling should occur between different configu¬ 
rations which satisfy the ice rules. This can be accomplished 
through the collective tunnelling of a group of protons, on any 
closed loop within the lattice, where the associated flux “ar¬ 
rows” also form a closed loop, as illustrated in Fig.|^ 

Where such a loop exists, it is possible to generate a sec¬ 
ond proton configuration satisfying the ice rules, simply by 
interchanging long and short proton bonds. This is equivalent 
to reversing the sense of all fluxes on the loop. The shortest 


loops for which this is possible in ice Ih consist of six oxygen- 
oxygen bonds, as illustrated in Fig. Quantum tunnelling of 
the form considered in this paper preserves the topological 
sector, since this is unchanged by any local rearrangements 
of protons which preserves the condition of local flux con¬ 
servation. As a consequence, the conserved flux is elevated 
to the role of a quantum number (441 . We note that, since 
all of these properties follow from the topological structure of 
ice, they are independent of the mechanism by which quantum 
tunnelling occurs 0001 . 

The effect of quantum tunnelling on loops of six bonds has 
previously been explored in the context of quantum dimer 
models on the diamond lattice (42^441 . and of quantum ef¬ 
fects in spin ice |[35l[37]|^|40l. These models have the same 
cubic symmetry, Fd3m, as Ic water ice, for which the shortest 
closed loop of bonds defines the edge of an hexagonal plaque- 
tte. Since all such plaquettes are related by lattice symmetries, 
a minimal model for quantum tunnelling can be obtained by 
introducing a single tunnelling matrix element g 

^tutelling = [I0)(0| + |0)(0|] , (3) 

o 

where the sum on o runs over all hexagonal plaquettes in the 
lattice, and tunnelling occurs between loops with the opposite 
sense of flux. We emphasise operators in Eq. act only on 
those plaquettes on which the flux arrows form a closed loop, 
as higlighted in Fig. These plaquettes are those described 
as “proton-ordered rings” in, e.g.. Refs. Ca and Co). 

In ice Ih, in contrast, there are two inequivalent types of 
six-sided plaquette, illustrated in Fig. The first is com¬ 
posed entirely of centre-symmetric bonds [cf. Fig. [^, and 
comprises the buckled hexagonal plaquettes which make up 
the hexagonal-symmetry layers of the Ih structure. A crystal 
with N oxygen atoms contains N/2 such plaquettes, which 
we label as type I [cf. Fig. [^a)]. The second is composed 
of four centre-symmetric and two mirror-symmetric bonds, 
linking neighbouring layers of the lattice. There are 3N/2 
such plaquettes, which we label as type II [cf. Fig.j^b)]. 

Since there are two, inequivalent, types of six-sided plaque¬ 
tte in the Ih lattice, the minimal model for quantum effects in 
ice Ih comprises two distinct matrix elements 

^nn^ng = “51 E [I I + I «) (O l] 

Oei 

-92 E [|0)(0| + |0)(0|], (4) 

Oeii 

acting on the space of all possible proton configurations obey¬ 
ing the ice rules. Symmetry alone does not place any con¬ 
straints on the values of gi and g 2 , and a priori, these matrix 
elements can take on either sign. 

It is important to note that collective coherent quantum tun¬ 
nelling of the kind described by Eq. 0 is distinct from in¬ 
coherent single-proton tunnelling. The study of proton tun¬ 
nelling in ice has rather a long history. The motion of single 
protons plays an important role in theories of ice’s dielectric 
properties (661 and it was believed for some time that this mo¬ 
tion may proceed by single particle quantum tunnelling 0 . 
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FIG. 4. (Color online). Flux representation of proton configurations in hexagonal (Ih) water ice. (a) Proton configuration, including a 
hexagonal plaquette of Type II (cf. Fig.|^. (b) Equivalent proton configuration represented using arrows. The displacement of protons (white 
spheres) from the midpoint of each oxygen-oxygen bond can be mapped to an arrow on that bond. The ice rules require that two “in” arrows 
and two “out” arrows meet at each vertex of the lattice (oxygen atom). Where arrows form a closed loop, it is possible to tunnel between 
different proton configurations satisfying the ice rules, by reversing the sense of all arrows on that loop. The letters A, D, and colour 
coding, indicate the convention for labelling oxygen sublattices adopted in Section [in| 


However, it is now generally believed that tunnelling of sin¬ 
gle protons is not effective in ice Ih at low temperature and 
ambient pressure | 62 l. Collective tunnelling of protons, how¬ 
ever, may still play a role in ice Ih, as it does in many other 
H-bonded systems (TEl I68lf70ll . The experimental results of 
Bove et al ca indicate that collective tunnelling of protons 
on rings of six H-bonds can indeed occur in ice Ih. 

Various ab initio studies have considered the possibility of 
proton tunnelling in water ice. Much of this work has been 
motivated by interest in the successive transformations be¬ 
tween cubic, high-pressure phases of ice, and in particular 
by ice X, an extreme high-pressure phase where the Bernal- 
Fowler ice rules no longer apply CD. Benoit et al ca find 
that the transformation from ice VIII to ice VII (which do obey 
the ice rules) is driven by quantum tunnelling of individual 
protons. This becomes more favourable under pressure, as the 
oxygen-oxygen distance decreases. A later study by Lin et 
al ca reinforces this conclusion, and highlights the emer¬ 
gence of collective, quantum tunnelling of groups of protons 
in cubic ice VII. Lin et al note that this type of tunnelling 
might also be effective in hexagonal ice Ih, albeit with a much 
smaller tunnelling matrix element. 

The case considered in this Article — collective, quan¬ 
tum tunnelling of protons in ice Ih — was recently studied 
in detail by Drechsel-Grau and Marx, using a combination 
of density-functional, molecular-dynamics and path-integral 
techniques ll^ . Considering a cluster of 48 water molecules, 
containing a single “proton-ordered” ring of the type shown 
in Fig.|^ Drechsel-Grau and Marx find that proton dynamics 


at low temperatures are dominated by the collective quantum 
tunnelling of protons from one ordered state of the ring to the 
other. A further study by the same authors El finds that par¬ 
tial deuteration of ice Ih will suppress this collective quantum 
tunnelling, in agreement with experimental results of Bove et 
al ca. Another recent ab initio study of nuclear quantum 
effects on the dielectric constant of water ice also finds strong 
quantum fluctuations of the proton system, consistent with the 
flndings of Drechsel-Grau and Marx El- 

If tunnelling of the form of Eq. 0 were to occur on a sin¬ 
gle isolated plaquette, the resulting ground state would be a 
quantum superposition of the two opposite senses of proton¬ 
ordering on the ring, in direct analogy with the resonating 
ground state of single benzene ring. However in a typical state 
obeying the ice rules, approximately a quarter of the plaque- 
ttes are “ordered”, and quantum tunnelling on these plaquettes 
allows the system to explore an exponentially large number of 
different states obeying the ice rules. 

In related three-dimensional quantum dimer and quantum 
spin ice models, dynamics of this type have been shown to 
stabilise a quantum liquid state, formed by a coherent super¬ 
position of an exponentially large number of states (351 Ell 
I39ll40ll42ti44ll . In water ice, such a state would be a massively 
entangled superposition of proton conflgurations obeying the 
ice rules, in which the molecular character of water molecules 
was preserved, but individual protons could no longer be as¬ 
signed to a given water molecule. This type of liquid should 
be contrasted with the “plaquette-ordered” phase found in 
simplifled two-dimensional models of water ice, which are 
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dominated by local resonances iSlCSHzll. 

At this time, we are not aware of any attempt to indepen¬ 
dently determine the two different matrix elements gi and g 2 , 

which define our model T^tunneNing initio 

simulation of ice Ih. And to the best of knowledge, the only 
available estimate of the magnitude of collective tunnelling in 
ice Ih is 


g - 1.46 X 10“"^ eV (5) 

taken from the density-functional calculations of Ihm HU. 
Intriguingly, this value is consistent with the energy scale of 
proton dynamics observed in inelastic neutron scattering d). 
And it is somewhat larger than the corresponding estimates 
for the putative “quantum spin ice” materials ll37]| . Taken at 
face value, this would suggest that water ice is potentially a 
more favorable place to look for the formation of an exotic 
quantum liquid state than are the quantum spin ices. 

With this in mind, in Section]^ of this article, we develop 
a theory of disordered proton configurations in the presence 
of quantum tunnelling, based on the minimal quantum model 
of ice Ih, ^tunnelling make the assumption that 

> 0 , ^2 > 0 , ( 6 ) 

so that the model is accessible to quantum Monte Carlo sim¬ 
ulation. Before examining the quantum model, however, it 
is necessary to understand the classical proton correlations 
which arise simply from the ice-rule constraint. This will 
form the subject of Section |I^ below. 


III. PROTON CORRELATIONS IN A CLASSICAL MODEL 
OF ICEIH 


In what follows, we develop a theory of proton correla¬ 
tions in a classical model of ice Ih, neglecting all quantum 
tunnelling between different proton configurations. This the¬ 
ory provides a detailed and microscopically-derivable phe¬ 
nomenology to explain the diffuse scattering which is arises 
as a result of static proton disorder in ice Ih (481SH EU EOl- 
[87l. 


We start, in Section [nrA| by developing a long-wavelength, 
classical field theory of proton configurations in ice Ih. This 
field theory has parallels with those developed to explain 
pinch-points in proton-bonded ferroelectrics C3 and spin- 
ice (541 . but displays a number of new features, which will 
become important in the quantum case. The details of these 
calculations are described in Appendix [B] 

In Section III B we introduce a lattice theory of proton cor¬ 
relations in classical ice Ih, and show that this reproduces the 
predictions of Section [HTA The details of these calculations, 
which are based on a generalisation of a method introduced 
for spin-ice by Henley El, are developed in Appendix 
Then, in Section [nrc| we use the lattice theory introduced 
in Section IIIB to make explicit predictions for the struc¬ 
ture factors measured in X-ray (80l [STl and neutron scatter¬ 
ing (48ll49l[8^ experiments. These show a number of inter¬ 
esting features, which we interpret in terms of the classical 


field theory developed in Section III A 


A. Continuum field-theory for protons in classical ice Ih 

The natural place to start in constructing a theory of pro¬ 
ton disorder in ice Ih is from the ice rules HlElia. The first 
ice rule states that each bond of the oxygen lattice contains 
exactly one proton [cf. Section |IJ . This proton is displaced, 
relative to the centre of the bond, towards one of the two oxy¬ 
gen atoms which make up the bond [cf. Fig[2. The displace¬ 
ment of this proton from the centre of the bond (rr') can be 
described using the Ising variable 

CFrr' = ±1 , (7) 

where arr' = +1 if the proton is displaced towards r', and 
(jrr' = —1 if it is displaced towards r. It follows that 

(Jrr' = —CFr'r- (8) 

The second ice rule states that each oxygen must form ex¬ 
actly two short (covalent) and two long (hydrogen) bonds with 
neighbouring protons. Written in terms of the Ising variable 
CFrr', this becomes a condition that, at every oxygen lattice 
site r 


TTlr — ^ ^ <^r+dri r — 0 (9) 

nn i 

where the sum runs over the four nearest—neigbours within 
the oxygen lattice, located at sites r + dri- For this purpose, 
it is necessary to divide the lattice of oxygen sites into four 
inequivalent sublattices as illustrated in Fig. These four sets 
of oxygen sites have different associated vectors dri, defined 
in Appendix [A| 

Just as in proton-bonded ferroelectrics ifT^ . or spin ice (541 . 
we can understand the proton correlations arising from Eq. 0 
by considering the spatial variation of the flux field repre¬ 
sented by the arrows on the bonds in Fig. For each 
oxygen-oxygen bond (rr'), we assign a flux 

Prr' = - - —CFrr' (10) 

ao 

where ao is the oxygen-oxygen bond distance. This fiux 
points in the direction of the polarisation of the associated H- 
bond, as illustrated in Fig. The total fiux from the four 
arrows around a single oxygen site is thus 

Pr = ^ ^ dr 2 CTp-l-dri r (H) 

— do 

nn i 

Knowledge of the fields rur [Eq. (|g] and Pr [Eq. ([g] on 
half of the oxygen sites (e.g. those on the A and C sublattices, 
shown in Eig. is sufficient to uniquely determine the proton 
configuration of the entire lattice. We can see this as follows: 
if, for a given oxygen site r, we know both rrir (which must 
be zero for an ice rule state) and Pr we can determine the 
value of all of the surrounding bond variables cTr+dri r using 
the relation 


<^r+dri r 



( 12 ) 
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Every bond belongs either to one A oxygen site or one C oxy¬ 
gen site [cf. Fig. [^b)] so knowing nir and Pr on just the A 
and C sites, (or equivalently just the B and D sites) is suffi¬ 
cient. 

We may imagine generating a proton configuration by set¬ 
ting mr = 0 on every A and C oxygen site and letting the 
fiux Pr vary between those sites. For the configuration thus 
obtained to be consistent with the ice rules we would need it 
to satisfy mr = 0 on all of the B and D tetrahedra as well. 
Naturally, this implies some constraints on the spatial varia¬ 
tion of Pr. These constraints control the form of the proton 
correlations. 

The constraints on the spatial variation Pr arising from 
the ice rules may be understood by defining continuum fields 
(i.e. defined over all space, not just on the lattice) Pa(i*) and 
Pc'(r) in such a way that evaluating them at the lattice sites 
returns the value of Pr [Eq. {TT}]. If we then assume 
Pa,c(i*) to vary smoothly in space, we can use the condition 
that mr must vanish at the B and D sites to obtain a constraint 
on the fields Pa,c{'^) and their derivatives. This procedure is 
described in more detail in Appendix]^ 

Neglecting terms beyond leading order in the bond distance 
ao, we obtain: 

-3Pl(r)+3P5(r)-2aoV-PA 

+3ao5,Pl - 3aod,P^ = 0 , (13) 

-3Pc{r) + 3P^(r) - 2aoV ■ Pc 

+3aod,PS - 3aod,PX = 0 , (14) 

where ^ is the hexagonal symmetry axis of the crystal (i.e. the 
crystalographic c-axis). These equations can be decoupled by 
introducing odd and even combinations of the fields Pa,c{'^) 

P+(r) = ^(P^(r)+P^(r)), (15) 

P_(r) = ^(P^(r)-P^(r)). (16) 

It follows that the uniform polarisation P+ (r) satisfies a zero- 
divergence condition 

V-P+(r) = 0, (17) 

while the staggered polarisation P-(r) is governed by the 
equation 

Pl{r) + ^aoV ■ P_(r) - aod,Pl{v) = 0 . (18) 

The very different form of the equations governing P+ and 
P_ suggest that these fields have qualitatively different cor¬ 
relations. Following (281 [54j [ 88 l, we can estimate these by 
assuming a free-energy of the form 

^ = (19) 

where k,^ is an (unknown) constant of entropic origin, and 
Vu.c. is the volume of a unit cell. Within this approximation. 


correlations of P± are controlled by a Gaussian distribution 
of fields 

d^T . 


p[P^] oc exp 


Ky f Pr 

“T J ^ 


( 20 ) 


subject to the constraints Eq. (Tf) and Eq. (18). 

We find that the correlations of P+ have the form of a sin¬ 
gular “pinch-point” in reciprocal space 

(-q)P^(q)) = T , 

{aA = x,y,z). (21) 

Meanwhile, the correlations of P_ in reciprocal space show 
much broader, smoothly varying structure 


1 




{a, 13 = x,y) 


where 


q_L = iqx,qy) 


( 22 ) 


(23) 


(24) 


(25) 


and, for compactness, we have introduced the notation 
^ ~ 3 ' 

Fourier-transforming Eq. ([2T|)-Eq. ( [2^ , we find that cor¬ 
relations of P+ decay algebraically in real space, with the 
dipolar form 


<5(r 




(26) 


are very short-ranged, decay- 
It follows that proton 


(P^(r)P^(O)) = ^ 

Meanwhile correlations of P_ 
ing over a length-scale C [Efi- ( [25] )] 
correlations at large distances are controlled by the field P+. 

The algebraic correlations of P+ give rise to sharp pinch- 
point singularities in structure factors, of the type observed 
by Li et al. in neutron scattering from ice Ih (48l . However 
in some Brilluouin zones, pinch-point singularities are sup¬ 
pressed by the lattice form-factor, and scattering is instead 
dominated by broad, assymetric features coming from the cor¬ 
relations of P_. We discuss this point further in Section III 
C, below, where we develop an explicit theory for neutron and 
X-ray scattering experiments. 

The form of the constraint on P+(r) [Eq. strongly 
suggests an analogy with electromagnetism, where the zero- 
divergence condition on magnetic field 


V-B = 0, 


(27) 


can be resolved as 


B = V X A , (28) 

and the electric and magnetic fields are connected by an under¬ 
lying 1/(1) gauge symmetry. This analogy, and the distinction 
between the two classical fields, P+ and P_, become explicit 
once quantum effects are taken into account, as described in 
SectionllVl 
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FIG. 5. (Color online). Correlation of long and short proton bonds in a classical model of hexagonal (Ih) water ice. The structure factor 
*S'ising(q) [Eq.|^, for the Ising variable a [Eqs. ([7])-(p^], is plotted in three orthogonal planes in reciprocal space. Near to zone centres, corre¬ 
lations are well-described by a combination of pinch-point singularities, reflecting the alg ebraic correlations of the held P +(q) [cf. Eq. j^], 
and smooth features, reflecting the short-range correlations of the held P-(q) [cf. Eq. ( |22|23| )], as discussed in Section [niA Calculations 
were performed using the method outlined in Appendix [C] in which the ice-rule constraints are written as orthogonality conditions in Eourier 
space El. Reciprocal-lattice vectors are indexed to the orthorhombic unit cell deflned in Appendix following the conventions of Nield 
and Whitworth 1491 • 


B. Lattice theory of proton correlations in classical ice Ih 


technical details for Appendix 


The classical fields P+ and P_, introduced in Section III A 


provide a complete description of the correlations of protons 
in classical ice Ih at long-wavelength, i.e. near to zone- 
centers in reciprocal space. However, X-ray and neutron scat¬ 
tering experiments on water ice measure proton correlations 
at all length-scales. We have therefore developed a lattice- 
based theory of proton correlations in classical ice Ih, valid 
for all wave-numbers. The approach we take is a generali¬ 
sation of the method developed for spin-ice by Henley m, 
in which the ice rules are expressed as a projection operator 
in reciprocal space. In what follows we explore how the pre¬ 
dictions of this theory relate to the those obtained from the 
continuum theory described in Section |III A| We reserve all 


The underlying structure of proton correlations in classical 
ice Ih is most easily understood throughout the correlations of 
the Ising variables cTyy' [Efi- which describe the alternat¬ 
ing long and short bonds between protons and neighbouring 
oxygen atoms. These are characterised by the structure factor 


'S'ising(q) = , (29) 

vv' 


where the sums on v and v' run over the eight distinct bonds 
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within the 4-site primitive unit cell, and 

CTi.(q) = \/ 4 X! ■ I^rr') 


rr'Gz^ 


— 


(30) 


In calculating crjy(q), we label the 4 oxygen sublattices within 
the unit cell A, B, C, D, (cf. Fig|^, and adopt a sign conven¬ 
tion such that 


rG{A,C}, T'e{B,D}. 


(31) 


In Fig.j^we show results for 5'ising(q), calculated within the 
lattice-based theory. The structure factor exhibits clear pinch- 
point singularities, characteristic of the ice rules OSIlll, at a 
subset of Brillouin-zone centers typified by 


q; = (o,o,2 ), 


(32) 


where, following Nield and Whitworth 1491 . we index all 
reciprocal-lattice vectors to the 8-site orthorhombic unit cell 
defined in Appendix [A[ 

Correlations near to reciprocal-lattice vectors (Brillouin- 
zone centres) are described by the classical field theory devel¬ 
oped in Section III A with contributions from both fields, P+ 


and P_. Near to a reciprocal lattice vector Q, for |q| ^ 1, 
the structure factor can be written 


^lsing(Q Q) 

« ^F;f'"s(Q)(|AQ;yp.(q)|2), 


(33) 


V = ± 


where the form-factor [Eq. ( |D12| )] and vectors 

^ Ising I I 

Aq ^ [Eq. (D13)] are define d in Appendix ID and the Fourier 
transform P^(q) [Eq. (C14)] in Appendix]^ 

Sharp pinch-points are seen for a subset of reciprocal lattice 
vectors Qp, for which 


Fi""S(Qp) = 0 . 


(34) 


In this case, correlations are controlled by the zero-divergence 
condition on P+(q) [Eq. (p^], and it follows from Eq. pT) 
and Eq. (|C14[) that 


5^lsing(Qp + Q) 


rising 


(Qp) 






Considering the specific example of 


for which 


This pinch-point, aligned with the hexagonal-symmetry axis 
(z-axis), can be clearly resolved in Fig.|^(b) and (c). 

While there are no reciprocal lattice vectors for which 
F^"^^(Q) vanishes identically with remaining fi¬ 

nite, there are another set of lattice vectors Qm, for which 




(39) 


and the structure factor is dominated by the short-ranged cor¬ 
relations of P_ [Eq. (18)]. Correlations of this type occurs 
for 


Q; = ( 1 , 1 , 0 ) 


(40) 


and are visible as a broad feature centred on this reciprocal 
lattice vector in Fig.[^(a). 

For a general reciprocal lattice vector Qpm 




(41) 


and correlations refiect a combination of pinch-points origi¬ 
nating in P + (q) and broad features originating in P_ (q). An 
example of this occurs for 

Q;m = (2,0,0) (42) 

visible in Fig. (a) and (c), where a pinch-point has been 
superimposed on a featureless background. 

Near to zone-centres, where they can be compared, the 
lattice-based theory is in complete agreement with the predic¬ 
tions of the classical field-theory developed in Section [III A[ 
In Appendix we show how the lattice-based theory reduces 
to the continuum theory at long-wavelength. We find the 
that entropic coefficient which controls correlations of the 
fields P^; [Eq. (^-^], is independent of i;, i.e. 


K, = = K,- 


(43) 


To confirm the validity of the lattice-based theory for more 
general q, we have also performed classical Monte Carlo sim¬ 
ulations of ice Ih, using local loop updates to sample proton 
configurations within the manifold of states obeying the ice 
rules. The results, described in Appendix]^ are in excellent 
agreement with the predictions of the lattice theory. 


C. Predictions for scattering from protons in a classical model 
of ice Ih 


X-ray and neutron scattering experiments on water ice 
do not measure the structure factor for bond-variables 

, but rather the 


5'ising(q) [Eq. (29)], discussed in Section 


IIIB 


Fourier transform of the correlation function for the density of 


we have 


q; = (0 ,0,2) 

(36) 

protons 




5'proton(q) = (n(-q)n(q)) . 

(44) 

Aq:®+ = (0 ,0,1) 

(37) 

Information about the proton disorder is contained in the dif¬ 
fuse part of this scattering, which is given by 18^ 



5'p)oton(q) = ®“(9l ■ ■ 

3-1/') , 

5ising(Q; + q)^ . 

(38) 

vv' 

(45) 
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FIG. 6. (Color online). Prediction for diffuse scattering of neutrons or X-rays from protons in a classical water ice described by the Bernal- 
Fowler “ice rules”. The structure factor Splotoni^) [Eq. (451], is plotted in three orthogonal planes in reciprocal space. The ice rules manifest 
themselves as both “pinch-points” — singular features in scattering, visible at e.g. Q* ,_|+ = (0,0,4) in (c) — and broad, asymmetric features 
at zone centers, visible at e.g. Q* |_|+ = (2,0,3) in (c). Results are shown for the theory described in Section ; 


me 


and Appendix C 


reciprocal-lattice vectors indexed to the orthorhombic unit cell defined in Appendix]^ following the conventions 6t iMeld and Whitwo; 


where cFjy{q) is given in Eq. (30), and is a set of vectors, 
defined in Appendix [A| such that the displacement of a proton 
from the midpoint on any given bond u is 


E)/y - 


with 


0.15 ao . 


(46) 


(47) 


In Fig. we show the results for diffuse scattering from 
protons in a classical model of ice Ih. The structure factor 
5p^;oto^^(q) [Eq. (45)], was calculated using the lattice theory 
introduced in Section |IIIB| At small momentum transfers, 
the scattering is suppressed by the factors of sin(q • a^y) in 
Eq. ([ 43 , and as a result, there is essentially no scattering for 
|q| < 2. For larger wave number, scattering shows a mixture 
of broad and sharp features, centered on two different sets of 
reciprocal lattice vectors. An example of broad feature can be 


seen near to Q = (2,0,3) in Fig.|^(c). An example of a sharp 
feature — a pinch-point — can be seen near to Q = (4,0,0) 
in Fig.[^(a) and (c). 

We can relate both broad and sharp features to the corre¬ 
lations of the classical fields P+ and P_, introduced in Sec- 
Expanding the structure factor 5'p^Son (^) [^4- ^45)] 
; reciprocal lattice vector Q, for |q| <C 1, we find 


III A 


tion 


about the 


cd if fuse 

O, 


^proton 


(Q + q) 


y]Fr-(Q)(iA^y"-p„(q)|2), 


(48) 


v=±: 


where the form-factors 


[Eq. (|D12|)] and vectors 


Aq ^ [Eq. (D14)] are defin ed in Appendi?^ and the Fourier 
transform P^(q) [Eq. (C14)] in Appendixj^ 

Once again, there are a subset of reciprocal lattice vectors 
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Qp,H + for which 




(49) 

and correlations are controlled by the zero-divergence con¬ 
dition on P+ [Eq. ([TtJ]. It follows from from Eq. ( [21] ) and 
Eq. (|C14[) that 


CiJ>"(Qp.H++q) 


T-iproton 

K 


(Qp,H+) 


1 - 


I C proton ^2' 




Considering the specific example of 
q;,h+ = (4,0,0) 

for which 


Ar°" ,.=( 1 , 0 , 0 ), 


we have 


cdiffuse 

^proton 


Qp,H+’^ 


(Qp,H+ + q) 


-I_^ 


(50) 

(51) 

(52) 

(53) 


A pinch-point singularity of this form is clearly visible near 
Q = (4, 0, 0) in Eig.|^(a) and (c). 

Similarly, while there are no reciprocal lattice vectors for 
which vanishes identically while re¬ 

mains finite, there are another set of lattice centers Qm,H +5 for 
which 


i"r‘°"(Qm.H+) « ^’r'°"(Q,,,H+) (54) 

and scattering from protons reflects the short-ranged correla¬ 
tions of P_ [Eq. ®]. An example of this type scattering 
occurs for 


Q;,h+= (0,4,1). (55) 

A broad, asymmetric feature, centred on this reciprocal lattice 
vector, is clearly visible in Eig.[^(b). 

Eor a more general choice of zone centre, Qpm,H +5 

n"'°"(Qpm.H+) « ^’r'°"(Qpn..H+) , (56) 

and scattering refiects the correlations of both P+ and P_. 
An example of this type scattering occurs for 


Q;m,H+= (0,4,3). (57) 

A combination of pinch points and broad, assymetric features 
can be seen near to this reciprocal lattice vector in Eig.[^(b). 

We conclude this discussion with a brief word of 
caution — in some cases, the scattering from protons 
^p^oton(q) 145)] exhibits pinch-points at the same re¬ 
ciprocal lattice vectors as pinch-points in the structure factor 
5'ising(q) [Efi- (29)] — as can be seen by comparing Eig.j^and 
Eig.[^ However, in general 

Ff°‘°"(Q)=0 Fi"'"®(Q)=0, (58) 


and pinch-points in S'pliSton (q) do not, necessarily, occur at the 
same reciprocal lattice vectors as pinch-points in 5'ising(q)- 
We will present a detailed comparison of these results with 
diffuse neutron scattering experiments on ice Ih in Section 

[Y3 


IV. PROTON TUNNELLING AND EMERGENT PHOTONS 
IN A QUANTUM MODEL OF ICE IH 


The analogy between the ice-rules and electromagnetism, 
underpinning the classical analysis of Section |I^ becomes 
complete once quantum effects are taken into account. In 
what follows, we show that the minimal model for quantum 
effects in ice Ih, [Eq. (|^], leads to a compact, frus¬ 

trated quantum U(l) lattice-gauge theory, with precisely the 
form of electromagnetism on a lattice. We explore the new 
features of a proton liquid described by such a theory, and 
make explicit predictions for both inelastic and quasi-elastic 
(energy-integrated) scattering of X-rays or neutrons from dis¬ 
ordered protons in a quantum ice Ih. This discussion proceeds 
as follows : 

Eirstly, in Section |IV A| we outline the derivation of this 
lattice-gauge theory. Technical details of these calculations 
are provided in App endix |E| 

Then, in Section |IV Bj we explore some of the features 
of lattice gauge theory in its deconfined (proton-disordered) 
phase. In particular, we show how its low-energy excita¬ 
tions can be thought of as the linearly—dispersing, birefrin- 
gent “photons”, and how these relate to the classical fields 
P+ and P_ introduced in Section |nl| 

Einally, in Section IV C| we discuss the experimental signa¬ 
tures of quantum water ice, described by the deconfined phase 
of the lattice gauge theory. 


A. Lattice-gauge theory 


Our route to a lattice-gauge theory of ice Ih closely paral¬ 
lels the cubic-symmetry case previously considered by Her- 
mele et al. 1351, and Benton et al Oil. The theory itself, 
however, contains a number of new features. 

We begin introducing a set of pseudospin-1/2 operators 
5'^, 5'+, S'” defined on the bonds rr' of the oxygen lattice. 
The ^(-component of the pseudo-spin is directly proportional 
to the Ising variable cJyv' [Eq. (0], introduced to describe 
proton-correlations in the classical case 


qZ _ _ _ _ qZ 


(59) 


In keeping with the directedness of [Eq. (59)] the ladder 
operators obey the identity 

5+, = 5^-;^ (60) 

The minimal quantum model for ice Ih [Eq. (|^], 

can be expressed in terms of these operators as 

<3i"ng = -9iY.[sts^sts^sts^+h.c] 


Del 

-92 Y. [StS^S+S^S+S^ 


h.c.] . (61) 


Oeii 


A mapping to a f/(l) lattice-gauge theory is then possible by 
writing the spin-1/2 operators in a quantum rotor representa¬ 
tion I 






S±, ^ 


(62) 
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subject to the canonical commutation relation 


The commutation relation Eq. ( [6^ is familiar in quantum 
electromagnetism as the commutation between an electric 
field E and a vector potential A, and substituting the rotor 
representation Eq. ( [6^ into Eq. (61) results in a compact U{1) 
gauge theory 

^compact ^ ^ ^ 

Oei 

-2g2 cos([Vo X A]) . (64) 

Oeii 

where the sum X]^rr')ecs centre-symmetric 

oxygen-oxygen bonds, X](rr')eMS mirror- 

symmetric bonds [cf. Eig. and [Vo x A] represents the 
lattice curl of A^r' around a hexagonal plaquette, which may 
be of type-I or type-II [cf. Eig.[^. 

The electric field E^r' subject to the constraint that 


Err' = ± - V bonds rr' . 
2 


(65) 


Eollowing ESI Ell, one may then argue that averaging 
over fast fiuctuations of the gauge field softens the constraint 
Eq. ( [65] ), and leads to a non-compact gauge theory on the links 
of the ice Ih lattice 

^2 E -®rr' + y E -^rr' 

(rr'>eCS (rr'>eMS 

+T E + T E ^66) 


Oei 


Oeii 


The parameters U and U' may be thought of as Lagrange mul¬ 
tipliers fixing the average value of E^^, on the two inequiva¬ 
lent types of bond. The average over fast fiuctuations will in 
general renormalise JC and 1C' from their “bare” values 


ICo — 2gi IC'q — 2^2- 


(67) 


On general grounds Esma, and by analogy with quantum 
spin ice ||35l[36l[39l, we anticipate that this lattice-gauge the¬ 
ory will possess both a deconfined phase, in which the protons 
form a disordered quantum fluid, and confined phase(s), in 
which the protons order. In what follows we confine our dis¬ 
cussion to the deconfined phase, without attempting to char¬ 
acterise any competing fixed points. 

Even with this restricted goal, the validity of 
7^u(i) [Eq. (66)] depends critically on the assumptions 
made in passing to a non-compact gauge theory. While 
these assumptions are reasonable, they can ultimately only be 
validated through quantum Monte Carlo simulation of the mi- 
croscopic model [Eq. (|6l](] — cf. OlllMlIlll and 

l44ll . We have therefore used variational quantum Monte 
Carlo (VMC) simulation of 'Hu(i) to establish that the de¬ 
confined phase of the lattice gauge theory, Hu(i) [Eq. ([66])], 
closely describes the correlations of the microscopic model 


^tunneNing [Efi- (61)], for the Symmetric choice of parameters 
gi = 92 - These results are presented in Appendix]^ 

In principle, it is also possible to extract the parameters of 
of the lattice gauge theory — U,U', 1C and 1C' — from de¬ 
tailed quantum Monte Carlo simulation of the microscopic 
model nulling [Eq- (61)] as a function of gi and g 2 — cf. 
Ref. ( 33 . However since the purpose of this Article is to ex¬ 
plore the properties of the deconfined phase, and no reliable 
estimates are yet available for gi and g 2 in real water ice, we 
will continue to treat U, U', 1C and 1C' as phenomenological 
parameters. 


B. Phenomenology of the deconfined phase : Why these 
photons ? 


In the absence of charges, the defining characteristic 
of the deconfined phase of the U{1) lattice gauge theory, 
Hu(i) [Eq. (66)], is its “photon”, a transverse excitation of the 
gauge field Awith definite polarization and linear dispersion 
at long wavelength. Since 'Hu(i) is quadratic in A, it can be 
solved by introducing a suitable basis for transverse fiuctua¬ 
tions of the gauge-field. This calculation is explained in detail 
in Appendix]^ following the methods described in Ref. (371. 
Here we concentrate instead on using this solution of 'Hu(i) to 
describe the new features which arise from the tunnelling of 
protons in water ice. 

In Eig. [^ we show the dispersion of the excitations 
of 7^u(i) [Eq. (66)], as they would appear in an inelas¬ 
tic neutron scattering experiment on ice Ih. The dynam¬ 
ical structure fa ctor f or coherent scattering from protons, 
5'coh(q, [Eq. (E21)], was calculated for the symmetric 
choice of parameters 


U=U' , 1C = 1C' , (68) 


and the dispersion has been normalised to the characteristic 
energy-scale of the lattice-gauge theory, \/UlC. Within this 
normalisation, the excitations have an overall bandwidth 


Acj _ fU' 1C'\ 
where, for this parameter set 


(69) 


/(1,1)^ 4.56. 


On the basis of published simulations for quantum spin- 
ice |[36H38]| . it is reasonable to expect that the bandwidth of 
the excitations of the gauge theory, Acc, should be of the same 
order of magnitude as the quantum tunnelling g. 

The excitations shown in Eig. [7] possess a number of strik¬ 
ing features, specific to ice Ih. At low energies, the model 
supports two, linearly-dispersing modes, with intensity which 
vanishes linearly approaching zero energy. These are the 
emergent “photons” of the lattice-gauge theory, and the fact 
that there are two such modes reflects the two possible polar¬ 
isations of the photon. The vanishing intensity of the photons 
at low energy is a feature shared with the emergent photons 
of (cubic) quantum spin ice, and refiects the fact that neutrons 
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FIG. 7. (Color online). Prediction for inelastic scattering f rom d isordered protons in a quantum model of ice Ih at zero temperature (T = 0). 
Results for the dynamical structure factor (q,^) [Eq. (E24i] are shown in three, orthogonal, planes in reciprocal space. The “photons” 
of the lattice gauge theory 77u(i) [Eq. ( [66| ] are visible as gapless, linearly-dispersing excitations at long-wavelength. These photons are 
birefringent, with a dispersion which depends on the polarisation of the photon, except along the optical axis, z. As a result two bands of 
photons are visible on the path F-M-K-T, while only one appears on the path F-A-L-H-A. Additional spectral weight at higher energies 
is associated with gapped, expo nentially-correlated excitations of the gauge-field. The dynamical structure factor for coherent scattering from 
protons (q, cj) [Eq. (E24i] was calculated using the U(l) lattice gauge theory 77u(i) [Eq. (66i], Results are shown for the parameter 
set given in Eq. (68), and were convoluted with a Gaussian of FWHM = O.OSbVWC to mimic the effect of experimental resolution. The 
color scale shows the intensity of the modes as they would appear in an inelastic scattering experiment, in the Brillouin zone centred on 
= (0,2,4). 


scatter from fluctuations of the proton density, and not directly 
from the underlying gauge fleld ifTTll . However, in contrast 
with the cubic case, the emergent photons of quantum ice Ih 
are birefringent, i.e. they have a dispersion which depends 
on the polarisation of the photon. The splitting of these two 
modes is clearly visible in Fig.|^ except for wavevectors par¬ 
allel to the hexagonal symmetry axis of the crystal, where they 
are degenerate. 

A second feature of note is the presence of gapped, optical 
modes. These are clearly visible in Fig. [7] at the zone center, 
r, at energies above the photon dispersion. There two such 
modes, and away from high symmetry points, they are gen¬ 
erally non-degenerate. The presence of these optical modes 
distinguishes quantum ice Ih from quantum spin ice, where 


the pure gauge theory only supports photons lJ7l . 


We wish to emphasize that the symmetric choice of parame¬ 
ters, Eq. ( [^ , is made purely for illustrative purposes, and that 
— as long as the lattice-gauge theory remains in its decon- 
flned phase — a more general choice of parameters will lead 
to qualitatively the same behaviour. In fact, as we shall show, 
the signature features of quantum ice Ih — the birefringence 
of the emergent photons, and the presence of gapped optical 
modes — are strongly constrained by symmetry, and follow 
naturally from the quantisation of the two classical flelds in¬ 
troduced in Section [Tin P+ P-* 


The correspondence between these classical flelds, and the 
excitations of the lattice-gauge theory Hu(i) [Eq. (66)], is a 
consequence of the fact that both theories automatically re- 
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FIG. 8. (Color online). Prediction for coherent, quasi-elastic scattering from protons in a quantum model of ice Ih at zero temperature 
(T = 0). The “pinch-points” associated with the ice-rules [cf. Fig. [^, are eliminated by quantum fluctuations [cf. Ref. 1^ 1. At finite 
temperatures, these pinch-points will be restored with a weight linear in T |32l. Results are shown for the energy-integrated, equal-time 
structure factor (q, t = 0) [Eq. 79 , calculated within the lattice-gauge theory 77u(i) [Eq. j66|)], for the same parameter set [Eq. ( 681 ] as 


Fig. El Reciprocal-lattice vectors are indexed to the orthorhombic unit cell defined in Appendix | A| following the conventions of Niela and 
Whitworth 14^. 


spect the ice-rules. It follows that constraints on the classical 
fields, Eq. ^TJ) and Eq. ( p^ , remain valid in the presence 
of quantum tunnelling [Eq. (|^], and that the long- 

wavelength dynamics of the quantum model can be found by 
quantizing the fluctuations of the classical fields. 

Let us consider first the case of P+. The constraint, 
Eq. can be enforced by writing 

P+ = V X A' (70) 

where it is important to distinguish the course-grained field 
A' from the microscopic field Apr', entering into the lattice 
gauge theory Hu(i). The form of the Lagrangian describing 
fluctuations of A' is then the Maxwell Lagrangian, subject to 
the hexagonal symmetry of the lattice. Choosing the Coulomb 


gauge 

V • A' = 0 (71) 


this is given by 



X A')“(V X A'f , 
(72) 

where the tensors pai 3 and are diagonal in the crystal ba- 
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sis, and have the form 



(73) 


(74) 


As a consequence, the photons described by [Eq. (72)] 
are degenerate when propagating with momentum parall^o 
the crystallographic z-axis, with dispersion 


UJ = 



(75) 


However, they are non-degenerate when propagating in the 
plane perpendicular to the crystallographic z-axis, with dis¬ 
persion 


u}i = .—\k\, U2 = J—\k\. (76) 

\ Pz \ P± 

which depends on the polarisation of the photon. 

We now turn to the field P_. This has exponentially, rather 
than algebraically, decaying correlations in the classical limit 
([^], and no associated gauge symmetry. 


[Eq. ( [22] ) and Eq. ___ 

However, because of the form of the constraint Eq. ( [T^ , fiuc- 
tuations of P_ can described by a Lagrangian written in terms 
of its planar component P^ : 


Cp_ = \jdtj - A(P^ 


±\2 


■ E E 

fiv'=x,y,z a^=x,y 




, (77) 


with the ^-component of P_ being fixed by Eq. (18). Once 
again, the form of the tensor is dictated by the symme¬ 
try of the lattice. Eor a general choice of we find that 
Lagrangian £p_ supports two gapped modes, which become 
degenerate approaching q = 0 


(78) 


In the light of this analysis, at long wavelength, we 
can identify the gapless, birefringent emergent photons of 
7^u(i) [Eq. ( 66 )] with quantised fiuctuations of P+, and its 
gapped opticaTmodes with quantised fiuctuations of P_. This 



correspondence is demonstrated explicitly in Appendix]^ 


C. Predictions for diffuse, coherent inelastic neutron 
scattering at low temperatures 

The formation of the quantum liquid state at T = 0 would 
have profound consequences for scattering experiments. In 
this section we will briefiy comment on what one could expect 
to observe in a measurement of the coherent scattering from 


water ice, in the scenario where quantum tunnelling of protons 
leads to the formation of a fluctuating U{1) liquid state. We 
shall postpone a consideration of the consequences for mea¬ 
surements of incoherent scattering, such as that performed by 
Bove et al GSl until our discussion of existing experiments in 
Section [Vj 

Eig. shows the equal time (energy-integrated) structure 
factor for the diffuse coherent scattering from disordered pro¬ 
tons at T = 0 

(q> i = 0) = 4 y] sin(q • ai) sin(q • a^) 

= 0)S|(q,t = 0)). (79) 

The pinch points are absent at T = 0, replaced by suppres¬ 
sions of the scattering around Brillion zone centers MM 

123 . 

This effect is most clearly understood by comparing Eig 
with the corresponding classical result shown in Eig. 
Around certain reciprocal lattice vectors, e.g. 

Qp = (2,0,0) (80) 

the classical scattering is directly proportional to a correlation 
function of the uniform polarisation P+, with no contribution 
from the staggered polarisation P_, and takes the form given 
in Eq. ( [53] ), i.e. 

C^o^nQp + q)^(i-^+^) . 

In the quantum case the pinch point form of Eq. ( [53] ) becomes 
modified by a factor of q, suppressing the pinch point 

5'h+(Q p + q, i = 0) oc g ^1 - j , (81) 

At finite temperature these pinch points are restored with a 
weight linear in tOU. This being the case, the clearest sig¬ 
nature of the formation of a quantum liquid which could be 
obtained from energy integrated scattering is the observation 
of a pinch point at high temperature, the intensity of which 
reduces in as the system cooled, heading towards a linear sup¬ 
pression of the scattering of the form of Eq. as T ^ 0. 
The nodal lines which are predicted in the classical scattering 
[Eig.|^ remain nodal in the quantum case [Eig.[^. 

Around other reciprocal lattice vectors, e.g 

Qm = (2,0,3) (82) 

there is a large contribution to the scattering from the fiuctu¬ 
ations of P^. In these cases the broad, assymetric features 
present in the classical scattering remain present in the quan¬ 
tum case at T = 0, but are now shifted to finite energy in 
accordance with the gapped nature of the fiuctuations of P^. 
Around reciprocal lattice vectors such as 

Qp^ = (0,4,3) (83) 

where the classical scattering shows a combination of pinch 
points and broad, assymetric features, the quantum theory pre¬ 
dicts that the pinch-point contribution will be linearly sup¬ 
pressed T = 0, as in Eq. ^T\ , while the broad feature will 
remain, albeit shifted to a finite energy. 
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This separation of the fluctuations of P+ and P_ as a func¬ 
tion of energy would be clearly manifested in a measurement 
of the inelastic scattering. This is illustrated in Fig. [7] which 
shows a prediction for the inelastic scattering around the re¬ 
ciprocal lattice vector Qjp = (0, 2,4) at T = 0. 

The linearly dispersing photon modes are visible with in¬ 
tensity / (X q (X uj, vanishing as they approach cc = 0 at the 
zone center. The gapped modes have finite weight approach¬ 
ing the zone center. Observation of these modes in an inelas¬ 
tic scattering experiment would represent convincing evidence 
for the formation of a protonic quantum liquid in ice Ih. 


V. APPLICATION TO EXPERIMENT 


In this Article we have developed a comprehensive theory 
of the disordered proton correlations in hexagonal (Ih) wa¬ 
ter ice, considering both a classical model based on the “ice 
rules” [Section]^, and a quantum model allowing for coher¬ 
ent quantum tunnelling of protons [Section |IV|. We anticipate 
that the classical theory should accurately describe proton cor¬ 
relations at temperatures where quantum effects are unimpor¬ 
tant, while the quantum theory becomes of interest at low tem¬ 
peratures, i.e. temperatures comparable with the tunnelling 
matrix elements 51,2 in [Eq. (Hi- 

In what follows, we place these results in the context of 
published experiment, exploring three particular themes : ice- 
rule correlations, as revealed by diffuse neutron scattering at 
relatively high temperatures [Sec. V A| ; low-energy excita¬ 
tions of protons at low temperatures, as revealed by the in¬ 
coherent inelastic neutron-scattering experiments of Bove et 
al. IT 8 I [Sec. rVBl; and the thermodynamics of ice at low tem¬ 
peratures [Sec. ra- 


A. Diffuse, coherent neutron scattering in the classical regime 


In Section IIIC of this Article we developed a theory of 
diffuse neutron scattering from hexagonal water ice based en¬ 
tirely on the Bernal-Fowler ice rules EEia. We found that 
the ice rules give rise to both pinch-point singularities, aris¬ 
ing from the algebraic correlations of the uniform polarisa¬ 
tion, P+ [Eq. (p~5])], and broad, assymetric features, arising 
from the short-ranged correlations of the staggered polarisa¬ 
tion, P_ [Eq. (p^]. The way in which these two fields con¬ 
tribute is controlled by form factors, and depends on the Bril- 
louin zone in question — cf. Appendix We now consider 
how these predictions, summarised in Fig. compare with 
experiment. 

We consider first the pinch points, originating in the uni¬ 
form polarisation, P+. That the ice rules give rise to pinch- 
point singularities in the structure factor is a very general 
result and widely known from the study of ice-like systems 
dESUMl. In the context of water ice, the prediction that 
the proton correlation function should be singular at Brillouin 
zone centres going back as far as Villain in 1972 EH. These 
pinch-point singularities are also visible in Monte Carlo sim¬ 
ulations of the structure factor, based on the ice rules EH 


and in reverse Monte Carlo fits to neutron scattering data 
ll4^l4^[83U85l . Other theoretical studies, which have utilised 
random walk approximations ESISSl and graph series expan¬ 
sions ISTII . while not specifically identifying the pinch points, 
have noted the presence of nodal lines in the structure factor. 

Experimental observation of pinch points in water ice is 
challenging, since coherent scattering of both neutrons and 
X-rays from protons is very weak, and the pinch points are 
located at reciprocal lattice vectors, where scattering is dom¬ 
inated by Bragg peaks associated with the ordered lattice of 
oxygen atoms. Nevertheless, pinch-point structures are visi¬ 
ble in the neutron-scattering from deuterated water ice, D 2 O, 
where coherent scattering is much stronger ll48l 1^ [Ml [84l . 

Clear examples of pinch points, occuring at the zone- 
centers predicted by our theory, can be seen in e.g. Fig. 2(a) of 
Li et al. ll48]| . or equivalently. Fig. 5(a) of Wehinger et al. ED. 
These should be compared with Fig. | 6 (c) of this Article, not- 
ing that these authors have indexed their reciprocal lattice vec¬ 
tors to an hexagonal unit cell, such that: 


H H L (hexagonal) 

^ 2H 0 L (orthorhombic) 


(84) 


The pinch-points present in the data are well reproduced by 
the theory developed in Section [TnCl Nodal lines, connecting 
the zone-centers where there are pinch points, can also be 
seen as marked suppressions of the scattering along certain 
high-symmetry directions in both theory and experiment. 

We now turn to the broad, assymetric features, originating 
in the staggered polarisation, P_. As well as pinch points, 
the neutron-scattering data of Li et al. also shows zone- 
center scattering with a broad, assymetric character. An ex¬ 
ample of a broad feature, arising from the correlations of P_ 
can be seen near to Qli = (1,1,5) in Lig. 2(a) of Ref.l48l or 
equivalently, Lig. 5(b) of Ref. [Ml This sho uld b e compared 
the scattering around Q = (2,0,5) in Lig. 6 (c) The broad 
zone-center features present in the data are well reproduced 
by theory. 

More generally, correlations at zone centers are described 
by a combination of pinch-points, originating in P+, and 
broad features, originating in P_. An example of a this type 
of feature, can be seen near to Qli = (2,0,3) in Lig. 3(a) 
of Ref. |48l or equivalently, Lig. 5(c) of Ref. [ST] This should 


be compared the scattering around Q = (0,4,3) in Lig. 6 (b) 
noting that 


77 0 1 / (hexagonal) 

0 2K L (orthorhombic) 


(85) 


Once again, there is good agreement between theory and ex¬ 
periment. 

Taking a broader view of reciprocal space, the proton corre¬ 
lations measured by neutron scattering from D 2 O at 20 K are 
generally well-described by the ice-rules, differing only in a 
stronger diffuse background, and “streaks” of diffuse scatter¬ 
ing in certain Brillouin zones where static proton correlations 
should not be visible Enillll- Both effects can be explained 
in terms of the thermal excitation of phonons ED- 

X-ray diffraction from protons does not suffer from the 
same problems as neutron scattering, and can be performed 
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FIG. 9. (Color online). Prediction for inelastic, incoherent scattering from protons in quantum ice Ih. (a) Prediction of the lattice-gauge 
theory, 7^u(i) [Eq. ([^], at a temperature of T = 5 K, for the same parameters set used in discussing coherent inelastic scattering [cf. Fig.l^ 
and Fig.|^. (b) Prediction of the lattice-gauge theory, convoluted with a Gaussian of FHWM 0.07 meV to represent finite experiments 
resolution, (c) Prediction of the lattice-gauge theory (blue), combined with the elastic contribution to the incoherent scattering (yellow). Both 
have been convoluted with a Gaussian to represent finite experimental resolution. The combined line-shape shows inelastic “wings”, similar 
to those observed in experiments on ice Ih by Bove et al d. Parameters are given in Eq. ( [68] ) and Eq. ( fShf , with details of calculations in 
Appendix [G| 


directly on H 2 O. However X-ray measurements present prob¬ 
lems of cooling, and to date all experiments on ice Ih have 
been carried out relatively high temperatures. In this case, 
diffuse scattering reveals relatively little about the ice rules, 
being dominated by a pattern of line-Iike “streaks”, charac¬ 
teristic of thermally excited phorons IMI . 


B. Incoherent inelastic neutron scattering at low temperatures 

At present, the strongest experimental evidence in support 
of the collective quantum tunnelling of protons in ice Ih comes 
from the recent neutron scattering experiments of Bove et 
al. im. In these experiments, inelastic, incoherent scatter¬ 
ing from protons in ice Ih at T = 5 K, was observed for a 
range of energies Acc ^0.1 meV, outside the experimental 
width of the elastic line. The authors found that the momen¬ 
tum dependence of this incoherent signal was consistent with 
a “double-well” model in which there is one proton on each 
bond, tunnelling between two sites. Since moving a single 
proton within a state obeying the ice rules has a prohibitive 
energy cost, they interpreted their result in terms of correlated 
tunnelling of protons on hexagonalplaquettes, of exactly the 
type described by [Eq. Q] — cf. Fig. ^ This in- 

terpretation finds support in ab inito calculations for Ih water 
ice (1911201, and the energy-scale of the dynamics observed 
is broadly consistent with the single published estimate of the 
scale of quantum tunnelling, g ^ 0.1 meV nSl. 

Given this, it is interesting to compare the predictions of 
the theory of water ice with quantum tunnelling developed in 
Section]^ of this Article, with the neutron-scattering experi¬ 
ments of Bove et al ca. Since incoherent scattering probes 
only the local proton correlations, it is not capable of distin¬ 
guishing the signal long-wavelength features of the lattice- 
gauge theory, such as birefringent photons, or temperature- 


dependent pinch-points. Nonetheless, the results of the com¬ 
parison remain very intriguing. 

In Fig. we present the results of a calculation of the in¬ 
coherent scattering at finite temperature T = 5 K, within the 


lattice-gauge theory 77u(i) [Eq. (66)]. Calculations were car¬ 
ried out for the for the symmetric choice of parameters 

u = u' ,k: = k:' 

[cf. Eq. ([^]. In the absence of further constraints on param¬ 
eters, we set 


'Juic = 0.018 meV , 


( 86 ) 


to give an overall bandwidth of excitations [Eq. l |6^ ] 

Aw ~ 0.1 meV ~ 1 K , (87) 

consistent with ah inito estimates of g (191 . The details of the 
calculation are given in Appendix [G| 

We find that the quantum tunnelling of protons does indeed 
produce “wings” of inelastic scattering which extend appre¬ 
ciably beyond the experimental width of the elastic line, as 
observed by Bove et al CD. Eine-structure in the incoher¬ 
ent inelastic scattering [Eig. |9(a)| , coming from the details of 
the dispersion [Eig.[7|, is obscured by finite experimental res¬ 
olution [Eig.|9(b)| , leading to broad wings on the elastic line 
[cf. Eig. 9(c)| . The energy-width of these wings is controlled 
by the bandwidth of collective excitatio ns of protons. Within 
the lattice-gauge theory this is set by \/WC and, for the pa¬ 
rameters chosen, is a little less than that observed in experi¬ 
ment. 

The good, qualitative, agreement between theory and ex¬ 
periment is very encouraging. None the less, it is hard to draw 
a definitive conclusion on the nature of ice Ih from incoher¬ 
ent scattering alone. A much cleaner test would be a mea¬ 
surement of the dispersion of the emergent photons, and as¬ 
sociated gapped modes, as a function of wave vector. This 
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would require coherent inelastic scattering, which is rendered 
rather challenging by the fact that incoherent neutron scatter¬ 
ing cross section for protons is approximately 50 times greater 
than the coherent cross section (891. However, in some cir¬ 
cumstances, it is possible to separate coherent and incoherent 
scattering using polarisation analysis (^ . 

It is possible to enhance the ratio of coherent to incoher¬ 
ent scattering by using deuterated ice, D 2 O (481 . However, 
some caution is called for. One of the key findings of the 
experiments of Bove et al. mil was that partial deutera- 
tion suppressed proton dynamics. This conclusion is sup¬ 
ported by subsequent ab initio simulations, which find that 
partial deuteration inhibits collective quantum tunnelling on 
“ordered” plaquettes El. 

It follows from the structure of the lattice gauge theory 
^u(i) [Eq. [^, that the quantum liquid is stable against all 
small perturbations which do not violate the ice rules (351. 
Consequently, the loss of quantum tunnelling on a very small 
proportion of plaquettes, through the natural abundance of 
deuterium in water, should not injure the quantum liquid state. 
Nonetheless, the loss of tunnelling on a macroscopic propor¬ 
tion of plaquettes is a different proposition, and could easily 
lead the protons to freeze. It is plausible that proton the dy¬ 
namics is restored in the case of full (or very large) deutera¬ 
tion. But even so the change from protons to deuterons will 
alter the relevant energy scale, presumably forcing quantum 
effects to lower temperatures, and changing the balance with 
any other interactions which favour proton order. 

A possible alternative to neutron scattering is X-ray diffrac¬ 
tion. In this it has been shown that at high temperatures the 
diffuse scattering is dominated by thermally excited phonons 
ED, making it difficult to observe the diffuse scattering 
which comes from the proton disorder. However, the quan¬ 
tum effects described in this article should manifest them¬ 
selves at temperatures far below the phonon Debye temper¬ 
ature. At these temperatures the contribution of thermally ex¬ 
cited phonons should be substantially reduced and one may 
hope to observe the correlations associated with the onset of a 
quantum proton liquid regime. 

It is also interesting to speculate that the optical properties 
of ice Ih, which so closely resemble those of the emergent 
photons in the lattice gauge theory, might be sensitive to a 
proton-liquid at low temperatures. 


C. Thermodynamics 

One of the most famous experimental results in the study of 
water ice is the demonstration by Giauque and Stout n that 
it retains a residual entropy down to T = I5K, and that the 
size of the residual entropy is very close to Pauling’s estimate 
of the entropy arising from disordered proton configurations 
obeying the Bernal-Fowler ice rules O. Subsequent experi¬ 
ments ||9TJ[92l have measured down to temperatures as low as 
T = 0.5K and find only a very small change in the entropy 
below I5K and no new features in the heat capacity. 

The quantum proton-liquid discussed in this article is a co¬ 
herent superposition of an exponentially large number of pro¬ 


ton configurations obeying the ice rules. At T = 0, it provides 
a unique quantum ground state, with vanishing entropy, in ac¬ 
cordance with the third law of thermodynamics (33. It fol¬ 
lows the residual entropy associated with ice must be lost in 
cooling from the high-temperature, classical, regime, to the 
low-temperature quantum regime. It is therefore natural to 
ask how this entropy would be lost, and whether this is con¬ 
sistent with published results for the heat capacity of ice Ih. 

In principle, entropy could be lost through either a 
sharp phase transition, or a smooth crossover. Estab¬ 
lishing which of these two possibilities occurs would re¬ 
quire finite-temperature quantum Monte Carlo simulations of 
^tunneNing [^9* which lie beyond the scope of the current 
Article. None me less, we can gain some insight by anal¬ 
ogy with quantum spin-ice, where related simulations have 
already been carried out. 

The nature of a thermal crossover between a classical and a 
quantum spin ice was explored in Ref. [571 The crossover was 
found to be controlled by a single length-scale Aj, the ther¬ 
mal de Broglie wavelength set by the thermal excitation of 
(emergent) photons. At small, but finite temperatures, the cor¬ 
relations return to their classical form at length scales r > Aj 
where 




( 88 ) 


The corollary of this result for scattering experiments is that 
the pinch point singularities, suppressed by quantum fiuctua- 
tions, are restored with a weight linear in Tl37l. 

The thermal excitation of photons also has consequences 
for the thermodynamics of a quantum ice. Since the photons 
of the lattice gauge theory are linearly dispersing, they give a 
contribution 


C - (89) 

to the heat capacity at low temperature (35] [571 EH- In the 
absence of a phase transition, this behaviour would be ex¬ 
pected to merge into a Schottky-like peak in the heat capacity, 
at a temperature of order the quantum tunnelling matrix ele¬ 
ment, g. 

The thermal-crossover scenario for quantum spin ice, de¬ 
veloped in Ref. [57] finds strong support in finite-temperature 
quantum Monte Carlo simulations (361 [381 • In particular, a 
recent numerical study by Kato and Onoda (551 presents a de¬ 
tailed analysis of the thermodynamics of a quantum spin-ice 
model as the temperature is lowered from a “classical ice” into 
a quantum spin-liquid state with emergent photon excitations. 
They found that the ice entropy is lost in a smooth crossover 
occuring at a temperature T < g, where g is the leading tun¬ 
nelling matrix element between ice states. This crossover is 
observable in the heat capacity as a low-temperature peak. At 
temperatures lower than the peak, the heat capacity behaves 
as C ~ T3, as expected for linearly-dispersing photons. 

Assuming that similar considerations apply, quantum ef¬ 
fects should begin to infiuence the thermodynamics of water 
ice at a temperature a little smaller than the tunnelling matrix 
elements gi^ 2 - Below this temperature, we would anticipate 
a peak in the heat capacity signalling the loss of the majority 
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of the Pauling ice-entropy. Below this peak there would be a 
substantial contribution to the specific heat from thermal- 
excitation of photons, in addition to usual the contribution 
coming from acoustic phonons. 

To date, to the best of our knowledge, measurements of ice 
Ih have not revealed any anomaly in the heat capacity which 
could be interpreted as the onset of quantum correlations of 
protons, down to T = 500 mK ll9^ . However, for a number 
of reasons, it is hard to draw any definitive conclusion from 
these experiments. 

Firstly, considerable ambiguity remains about the temper¬ 
ature at which quantum effects should be expected to occur 
in ice Ih. The single published ab initio estimate of quantum 
tunnelling ^ ^ 0.1 meV ^ 1 K |[T9ll . is very similar to the 
lowest temperatures achieved in experiment 19^ . Moreover, 
in the absence of finite-temperature quantum Monte Carlo 
simulations, it is not known at what fraction of gi and g 2 a 
heat-capacity anomaly should be expected to appear. It may 
therefore be that quantum tunnelling of protons does lead to 
quantum-liquid state in ice Ih, but that the associated loss of 
entropy occurs at a temperature lower than that currently mea¬ 
sured. 

The second problem associated with the interpretation of 
theremodynamic measurements low temperatures is the ex¬ 
treme difficulty of performing experiments, in equilibrium, on 
a system which retains an extensive entropy. Here, parallel 
studies of spin ice provide a stark warning : recent experi¬ 
ments on the spin ice Dy 2 Ti 207 by Pomaranski et al revealed 
thermal equilibriation times in excess of a week at tempera¬ 
tures of order 300 mK (931. Once their sample had equilib- 
riated, Pomaranski et al found an upturn in the heat capac¬ 
ity at low temperatures, in contrast with earlier experiments, 
where the Pauling ice-entropy was reported to persist down 
to 200 mK I8I9I . If similar problems of equillibriation occur, 
it is possible to envisage that a specific heat peak associated 
with the loss of Pauling’s ice-entropy could yet be observed 
in ice Ih at temperatures of order 1 K. 

While thermodynamic measurements at very low temper¬ 
atures will always be challenging, it is important to note that 
the consequences of quantum tunnelling in water ice should be 
observable at temperatures considerably greater than that as¬ 
sociated with the loss of the Pauling ice-entropy. This is par¬ 
ticularly true of dynamical properties, measured at wavelength 
shorter than the thermal de Broglie wavelength, Aj. As a re¬ 
sult, Quantum Monte Carlo simulation of quantum spin ices 
show clear signs of quantum effects in dynamical structure 
factors at relatively high temperatures |[36l[38l[78l. It therefore 
remains reasonable to discuss the inelastic neutron scattering 
of Bove et al CD, carried out at 5 K, in term s of a lattice 
gauge theory with characteristic energy scale s/UlC ~ 1 K 
[cf. Section|vg. 


VI. CONCLUSIONS 

Common, hexagonal, water ice is a wholly remarkable sub¬ 
stance — a proton-bonded network of water molecules in 
which oxygen atoms form a regular crystal, but protons need 


never order. In this Article we have explored the nature of pro¬ 
ton disorder in hexagonal (Ih) water ice, considering both a 
classical model based on the Bernal-Fowler “ice-rules” CH, 
and a quantum model which respects the ice rules, but also in¬ 
corporates collective quantum tunnelling of protons on hexag¬ 
onal plaquettes [cf. Fig.|^. Quantum tunnelling of this type 
is known to have profound consequences in models of quan¬ 
tum spin ice (35^401 . and there is a growing body of evidence 
that it also plays a role in water ice, including ab initio calcu¬ 
lations (131201, inelastic neutron scattering ca, and parallel 
studies of other proton-bonded systems (T^lbSl . 

In the case of classical ice Ih, we have developed a com¬ 
prehensive theory of diffuse scattering from protons [Sec¬ 
tion |I^. We find that the ice-rules have two distinct signa¬ 
tures in scattering : singular “pinch-points”, originating in 
a zero-divergence condition on the uniform polarisation P+, 
and broad, assymetric, zone-center features, coming from the 
staggered polarisation P_ [Fig. [^. Both of these features 
have previously been observed in experiment (481 . 

In the case of the quantum model, we have obtained 
a description of a low-temperature quantum liquid — in 
which protons resonate between an exponentially large num¬ 
ber of configurations satisfying the ice rules — in terms of a 
quantum U{1) lattice-gauge theory [Section]^. The long- 
wavelength excitations of this quantum liquid take the form 
of gapless, emergent photons, originating in the uniform po¬ 
larisation P+, and gapped, optical modes originating in the 
staggered polarisation P_. 

We have used this lattice-gauge theory to make concrete 
predictions for inelastic scattering experiments on a ice Ih 
[Fig. [71. We find that both the emergent photons and the op¬ 
tical modes can be clearly resolved at finite energy. Much 
like real light in water ice (94j[95l, the emergent photons are 
birefringent, exhibiting a dispersion which depends on both 
their polarisation and their direction of propagation. We have 
also explored how quantum tunnelling of protons modifies 
diffuse scattering at low temperatures [Fig.j^. We find that 
the “pinch-points”, characteristic of the ice-rules, are progres¬ 
sively eliminated as the system is cooled toward T = 0. 

The assertion, in Ref. EH, that the lattice gauge theory 
describing a two-dimensional quantum water ice should be 
confining at T = 0, is consistent with numerical results for 
two-dimensional quantum ice (TSl [7614781 . However a num¬ 
ber of three-dimensional ice-like models are known to sup¬ 
port deconfined, quantum-liquid ground states (33 SH. In the 
absence of detailed simulations it is not clear for which pa¬ 
rameters the ground state of 'Htunneiiing should be or¬ 

dered. However the incoherent inelastic neutron scattering ex¬ 
periments of Bove et al (T^ provide prima facie evidence of 
collective proton dynamics at a temperature (5 K) comparable 
with ab initio estimates of quantum tunnelling in ice Ih CD, 
and is in good, qualitative, agreement with the predictions of 
our theory [Fig.j^l. 

The observation of a birefringent, emergent, photon in co¬ 
herent inelastic scattering from protons at low temperatures 
would provide very strong evidence for the existence of quan¬ 
tum fiuid of protons in ice Ih. At present, however it is not 
possible to compare these predictions directly with expert- 
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ment, since coherent inelastic scattering data is unavailable 
for the temperatures where quantum tunnelling is expected to 
relevant. Experimental evidence for a quantum fluid of pro¬ 
tons from thermodynamic measurements is also lacking [cf. 
Section |V|. We hope that further experiment will provide a 
deflnitive answer to these questions. 

Experiments probing the excitation spectrum of other water 
ices, where protons order at low temperatures, could also be of 
interest, since experience with other ice models suggests that 
a quantum liquid may still be observable at flnite temperature 
(361 Ell UM • And in the absence of evidence to the contrary, 
a quantum-liquid ground state in ice Ih remains an intriguing 
possibility. 

In conclusion, more than 80 years after the pioneering work 
of Bernal and Eowler, the behaviour of protons in common, 
hexagonal water ice at low temperatures remains a problem of 
great fundamental interest. And, on the basis recent experi¬ 
ments ca, a quantum liquid of protons, of the type explored 
in this Article, remains a tantalising possibility. 
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Since this work was completed, two papers of direct rele¬ 
vance have appeared. The first, by Isakov et al (9611 , also 
examines the proton correlations in a classical model of ice 
Ih. The second, by Yen and Gao (971 , reports evidence for 
quantum coherence in water ice, from measurements of its 
dieletric constant. 


Appendix A: Details of the lattice and coordinate system 

In this Appendix we give details of the lattice and coordi¬ 
nate system which we have used in our calculations. We will 
write all lattice length scales in terms of the oxygen-oxygen 
bond distance, ao, which in ice Ih is approximately m 

ao « 2.75 A. (Al) 


The primitive unit cell contains four oxygen atoms and eight 
protons. 

One can also define an orthorhombic unit cell, containing 
eight oxygen atoms with orthogonal repeat vectors 


Gx = Gi (A5) 

Gy = 2 G 2 - Gi (A 6 ) 

Gz = G 3 . (A7) 


This is the unit cell used by Nield and Whitworth in Ref. (49l . 

In the scattering patterns shown in Pigs. [^|^[^and[^ the 
reciprocal lattice units used for the momentum scale are with 
reference to the orthorhombic unit cell, as in Ref. (4^ . Thus, 
h, k and I [cf. Pigs. |^[^[^and[^ relate to the momentum 
transfer q via 


q = 


27rh 2'Kk 27 tI 


'^1 


IG; 


(A8) 


where Gxyz are defined in Eqs. (A5)-(A7). 

The bond vectors df connecting an oxygen to it’s four 
neighbours [cf. Eq. ([^] are given by 

=d[( =ao (0,0,1) (A9) 

(AlO) 


= ao \ X, 


da = ao - 


dg = ao 



df = ao I 


(All) 


(A12) 


(A13) 


df = ao 



3 ’ 


d3"^=ao (^0,?|^, 

df = —df V i 
df = -df V. 


(A14) 

(A15) 

(A16) 

(A17) 


The displacement of protons from the bond midpoint on the 
8 sublattices of bonds are given by where 


We use a coordinate system in which the hexagonal symme¬ 
try axis is the z-axis and the repeat vectors of the hexagonal 
unit cell 


ao = 0do 

a4 = 0df 


ai = 4)A^ 3l2 = </>df as = 4)A^ 
as = (/)df ag = ^df sq = ^df . 

(A18) 


Gi — ao 



G2 = ao 



Gs 


ao 



(A2) 

(A3) 

(A4) 


Por the purpose of our scattering calculations we take the size 
of the proton displacement from the bond midpoint relative to 
the bond length to be 

0 = 0.15 (A19) 

assuming an 0-H covalent bond length (Ml 
<^OH = 0.95 A. 


(A20) 











(a) Bond labelling convention 


(b) Plaquettes labelled ‘0’ and ’4’ 


(c) Plaquettes labelled ‘1’ and ’5’ 



(d) Plaquettes labelled ‘2’ and ’6’ (e) Plaquettes labelled ‘3’ and ’7’ 

FIG. 10. (Color online). Labelling convention for the eight sets of oxygen-oxygen bonds and eight sets of six-link plaquettes unrelated by 
translational symmetry. This is the convention employed in diagonalising the lattice gauge theory described in Section[lV|and Appendix]^ (a) 
is color-coded such that bonds related by translational symmetry are the same color. The bonds labelled 0 and 4 are mirror symmetric bonds 
aligned along the z-axis and the remaining bonds are centre symmetric [cf. Fig. |^. (b) shows plaquettes of type I with matrix element gi. 
(c)-(e) shows plaquettes of type II with matrix element g 2 . 


For the calculation of the dispersion of the lattice gauge the¬ 
ory, presented in Section]^ and Appendix]^ it is necessary 
to define a labelling convention for the eight sets of oxygen- 
oxygen bonds not related by translational symmetries and for 
the eight sets six-link plaquettes. Our convention is defined in 
Fig.[^ 

Of particular are importance are the set of vectors which 
link the bond midpoints around the edge of a plaquette to the 
center of the plaquette. We denote these vectors c^m where 
two bonds of sublattice m belong to a plaquette of sublattice 
n and are located at ± c^m with the plaquette center being 
at Pn. Where only one bond of sublattice m belongs to a 
plaquette of sublattice n we denote it’s position relative to the 


centre of the plaquette by C^m- These vectors are 


coi = —cio = C45 = —C54 = ao 

C02 = —C20 = C46 = —C64 = CLO 

Cos = —Cso = C47 = —C74 = ao 


( 4 , 0 , 0 ) 


C12 — C21 — ao 



42 A 

3 ’3) 


Ci3 — C31 — ao 



1 2 
3 
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C23 = ^32 = CiO 
C 56 = Ces = ciQ 
C 57 = C 75 = ao 
C^67 = C 76 = ao 
C 16 = C 25 = ao 
Ci 7 = C 35 = ao 
C 27 = C 36 = ao 
C 53 = C 71 = ao 
C^63 = ^^72 = ao 
Cei = C 52 = ao 


Vd' 3 ^ 2 ’ 3 



\/6’ 3 V 2 ' 3, 

1 1 2 ^ 



3 V 2 ’ 3 

1 1 


3v^’ 3 

1 1 2 
\/6’3\/2’ 3 

1 1 2 
7^’ 7^’ 3 



Appendix B: Derivation of constraints within in continuum 
field-theory 


In this Appendix we describe in detail the coarse-graining 
procedure used in Section |III A| to obtain the constraint 
Eqs. ([^-([^. 

The argument begins by setting mr = 0 for all A and C 
vertices [cf. Fig. |^b)]. If we then set the value of the total 
flux Pr around every A and C vertex, we have specifled the 
H-bond conflguration of the entire lattice. The necessity that 
this conflguration must also obey mr = 0 at the B and D 
vertices induces constraints on how Pr varies in space. 

To derive those constraints in the long wavelength limit we 
introduce two coarse-grained flelds P"^(r), P^(r) deflned at 
every point in space (not just on the lattice). We deflne P"^(r) 
in such way that when it is evaluated at the position of an 
oxygen vertex of the A sublattice, it returns the precise value 
of the total flux around that vertex 


P^(r = r^) =Pr^. (Bl) 

Similarly 

P^(r = rc) = Prc- (B2) 

To obtain our constraints on P"^(r), P^(r) we write the 
scalar held mr at the B and D vertices in terms of the degrees 
of freedom at the four surrounding vertices. 

A vertex of the B sublattice [cf. Fig. [^, located at r^, 
neighbours one A vertex located at 


.( 0 ) _ 


= - do 


(B3) 

(B4) 


and three C vertices located at 

r^)=rB-df i = 1,2,3 (B5) 

[cf. Eqs. (|A^-( |AT71 i]. 

Using the vanishing of mr at the A and C vertices 

and similarly 

= ! {E } ■ <B7) 

Assuming smooth variation of P"^(r), P^(r) we can use a 
Taylor expansion to write 

Pr\c+Sr ^ Plcir) + . VP^c(r). (B8) 

Enforcing mr^ = rurj^ = 0 throughout the lattice, we ob¬ 
tain the constraints on our continuum fleld-theory, Eq. tH- 

•El- 


Appendix C: Projection operator calculation of structure factor 
for classical water ice 


Here we show how to calculate the structure factor for the 
Ising bond variables a^r' [Eq. ([29| )] via a generalisation of 
the lattice calculation in Ref. (5^ 7 In this approximation the 
Ising nature of the variables cTyv' is relaxed, such that they can 
take on any real value and their normalisation is enforced on 
average 


{al,) = 1 . 


(Cl) 


The calculation proceeds by constructing a projection opera¬ 
tor which acts on the Fourier transform of an arbitrary proton 
conflguration, to remove all states which do not satisfy the ice 
rules. 

We begin by deflning the Fourier transform over the bond 
variables as in Section ing [Eq. ([3T])], i.e. 


(q) = y ^ ^ 




H / — 


where v indexes one of the eight sublattices of bonds which 
are not related by a translational symmetry of the lattice, and 
N counts the number of oxygen atoms. The sign of cF^r' is 
flxed by the convention [Eq. ( [3T] )] 

rG{d4,C}, v' e{B,D}. 


where a = A^B^C^D indexes the four oxygen sublattices 

(cf. Fig. 13b)). 
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From these Fourier transforms we define an 8 -component 
vector 

o-(q) = 

(o-i (q), 0-2 (q), 0-3 (q), (74 (q), (75 (q), £76 (q), (77 (q), (78 (q)) . 

(C2) 

We also introduce a net Ising polarisation for each oxygen 
sublattice a = A, D, 

rriaiq) = ' r)m(r) . (C3) 

rGa 

It follows from Eq. ([^ that 

ma{q) = 0 Va,q (C4) 

in any state obeying the ice rules. 

The constraint Eq. can be written as a set of orthogo¬ 
nality conditions on d(q) 

(T'a(q),5-(q)) = 0, a = A,B,C,D (C5) 

where the inner product 

= (C 6 ) 

i 

and the 8-component vectors Tq, (q) are defined by 

ta = 


(giq.d„^/2^ giq.d^/2^ giq-d^/2^ giq.d3^/2^ g^ g^ 


TB = 


(giq-d®/2^ g^ g^ g^ g^ giq-df/2^ giq-df /2 ^ /2^ 


TC = 


( 0 , 0 , 0 , 0 , 


(C7) 


(C 8 ) 


(C9) 


td = 


( 0 , / 2 ^ 0 , 00 ) 


(CIO) 


with the vector df defined in Eqs. ( |A9| l-( AT7) |. 

At long wavelength (i.e. in the vicinity of a Brillouin zone 
center), the orthogonality conditions Eqs. (|C7li-(|CT0li must 
reduce to the constraints on the classical fields P+ and P_, 
discussed Section |III A| and Appendix The connection be¬ 
tween the two approaches can be made explicit by rewriting 
the constraints on bond variables, Eqs. ( |C7| )-( [CT0| ), in terms 
of the lattice variable nir and Pr. To this end we introduce 
the Fourier transform of the polarisation Pr on subblatice a 


Pa(q) = y ^ X]exp(-iq- r)Pr , 


(Cii) 


For the moment, we restrict our discussion to wavevectors 
q which are close to a reciprocal lattice vector Q and write 

q = Q + q (C 12 ) 

By analogy with Eqs. of the main text, we write 

m±{q) = A [exp(iQ • r^) mA{q) 


±exp(iQ ■ rc) toc(Q + q)] 


(C13) 


P±(q) = [exp(iQ ■ r^) Pa(Q + q) 


±exp(iQ • rc) Pc(Q + q)] • 

(C14) 

The vectors and vc are the positions of the oxygen sites 
labelled A and C within a primitive unit cell. Note that we 
have defined the fields in Eqs. ( [C13| ) and ( |C14| ) in such a way 
that they are independent of the reciprocal lattice vector Q 
and depend only on the distance to the zone center q. We may 


therefore use Eqs. (|^-(|23]), with the replacement 

q ^ q 


(C15) 


to describe the correlations at small q in the vicinity of all re¬ 
ciprocal lattice vectors Q, i.e. near all Brillouin zone centers. 

To demonstrate that the constraints derived in the con¬ 
tinuum theory [Eqs. are equivalent to the con¬ 

straints in the lattice theory [Eq. ( |C5| ) ] in the limit ^ ^ 0 
we consider the following linear combinations of the vectors 
r(q) = r(Q + q) 

TAC+(Q,q) = ;^(exp(iQ-rA)TA(Q + q) 

+(exp(zQ • rc)Tc(Q + q)) (C16) 

TAC-(Q,q) = ;^(exp(zQ • rA)TA{Q + q) 


- exp(iQ ■ rc)Tc(Q + q)) 


(C17) 


7 -BD+(Q,q) = ;^(exp(iQ-rB)TB(QH 
exp(iQ ■ rD)TD{Q + q)) 


TBD-iQA) = ;^(exp(iQ ■ rB)rB(Q - 
exp(iQ • rB)TB(Q + q)) 


q) 


(C18) 


q) - 


(C19) 


where and r u are the positions of each oxygen 

vertex in the unit cell. 


rGa 


We can then express the constraints on d(q) [Eq. (C5)] as 

(rAc+(Q,q),CT(Q + q)) = 0 (C20) 

(r^C-(Q,q),a(Q + q)) =0 (C21) 

( 7 BD+(Q,q),CT(Q + q)) = 0 (C22) 

(TBB-(Q,q),CT(Q + q)) = 0 . (C23) 
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Since 


(r^C+(Q, q), 5'(Q + q)) = TO+(q) (C24) 
(r_Ac-(Q,q),5'(Q + q)) = «^-(q), (C25) 


the constraints Eq. ( |C20| l and Eq. ( |C21| > are satisfied if 

m+(q) = m_(q) =0 V q . (C26) 

The constraints on r(Q, q), Eq. ( [C22| ) and Eq. ( [C23| ), give 
rise to constraints o n the q dependence of P+ (q) and P_ (q). 
Reexpressing Eq. \C22 ) in terms of P+(q) [Eq. ( [Tsj i] and 
P_(q)[Eq. @], and expanding to linear order in q we find 


q'P+(q) + c>(ao) = 0 


(C27) 


This precisely the Eourier transform of Eq. JTt]., i.e. the re¬ 
quired condition on P+. Similarly, Eq. ( C23| ) becomes 

PI (q) + iao Qq • P" (q) - q^Pl (q)^ + O(a^) = 0 

(C28) 


which is the Eourier transform of Eq. ( p^ . 

Having established that the lattice-based theory is equiva¬ 
lent to the continuum field theory at long wavelength, we now 
turn to the problem of calculating the structure factors which 
describe proton-proton correlations for arbitrary q. Within 
the lattice-based theory, this reduces to constructing a matrix 
V which projects states into the subspace of proton configura¬ 
tions orthogonal to the set of vectors (q). Explicit construc¬ 
tion of V is messy, and the final expression for the structure 
factor must be evaluated numerically. 

We proceed by using Gram-Schmidt orthogonalisation 1981 
to construct from {ra{q)} an orthogonal basis set {^^(q)}' 


't'a = TA 


Tb=tb- 

TC = TC- 


(t'a^'^a) 


t'a 


Ta - 


(tct'b) r 


rD=rD- 


, {td,Tb), 




Ta 


(t'e^Tb) ' 


(C29) 

(C30) 

(C31) 


{td,Tc) 


(C32) 


The projection matrix P(q) is then given by 


Pii(q) = Sij - 

a 


‘ Oil' OiJ 


(C33) 


The structure factor is obtained by acting on the Eourier 
transform of some general (non-ice rule obeying) state d'(q) 


(cri(-q)c7j(q)) = 

= X]Pm(-q)Pim(q)«(-q)o-;„(q)) 

mn 

~ •V)^mnPin{ q)Pjm(q) 

mn 

= MaViM) (C 34 ) 


where A/q is a normalisation constant and in the last step we 
have used the identities 

P(-q) = P(q)* (C35) 

P(q)* = P{^f (C36) 

P(q)2 = P(q). (C37) 

The constant A/q is fixed by Eq. This method was 

used to calculate the scattering patterns in Eigs. [^an dm 


Appendix D: Relating structure factors to correlations of P+ (q) 
and P-(q) 

In this Appendix we derive the relationship between the 
structure factors 5'ising(q) and discussed in Sec¬ 

tion [IIT^ and Section [HTcI of the main text, and the correla¬ 


tion functions of the fields P±, introduced in Section III A 


In so doing we establish the necessary conditions for a zone 
centre to exhibit singular, pinch-point scattering. 

The structure factors 5'ising(q) [Eq. (29)] and 
tt(q) [Eq. (|45l] can both be expressed in terms of 


the correlations oftne Ising variable crjy(q) [Eq. ([^], 


as 


5'x(q) = y](o-r^(q)o-r^'(-q)) vU^) (di) 


where 

X = Ising , proton , (D2) 

the sum on a runs over all bonds within the unit cell, and the 
coefficients r]^ (q) depend on which structure factor is being 
calculated. In the case of 5'ising(q) 

V'r^ici) = l Vi.,q, 

while for S'^^oton (q) 

.P-t-(q)=zsin(q.a,), 


(D3) 


% 


(D4) 


where the vectors describe the displacement of the protons 
from their bond midpoints, as defined in Eq. ( |A18| ). 

The expression for S^{q) [Eq. (Dl)] can be factorised to 
give 


'5x(q) = (| y]^J(q)o'i'(q)r 


(D5) 


We are interested in understanding the behaviour of this struc¬ 
ture factor near to a given reciprocal lattice vector Q, in terms 


of the continuum field theory developed in Section |III A| To 
this end, we write 


q = Q + q 


(D6) 


E]»7?(Q + q) o-i.(Q + q) 

V 

= Ai+(Q,q) rn+{q) + iP{Q,q) m_(q) 

+ A^(Q,q)-P+(q)+ A^(Q,q)-P_(q) , 

(D7) 
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whe re the vectors are defined in Eq. (D13) and 

Eq. (D14), below, and the scalar functions /i±(Q, q) play no 
part in our subsequent discussion. 

In any state obeying the ice rules 


m+(q) = m_(q) = 0 (D8) 


so terms in these fields can safely be dropped 

E]»7?(Q + q)o-i.(Q + q) 

V 

= A^(Q,q)-P+(q)+A^(Q,q)-P_(q) . 

(D9) 

It follows tha t, for |q| ^ 1, the structure factor 

S^{q) [Eq. |D5| can be written as 

SxiQ + q) « E • P-(q)l') ’ (DIO) 

i;=± 


[Eq. ([45])], we have 

Ar°"(Q,q) 


= exp(-zQ • r^) E exp[-z(Q + q) • df/2] 


i=0 


xisin[0 (Q + q • df )] 


±exp(-iQ • vc) ^ exp[-i(Q + q) • df/2] 


i=0 


xi sin[0 (Q + q • df)] 


(D14) 


where the parameter 

0 = 0.15 (D15) 

expresses the relative displacement of the protons from the 
midpoint of the bond, as defined in Eq. ( |A19| ). 


where Aq ^ is a unit vector in the direction of 


Aq,„ = AJ(Q,q = 0) , 

[cf. Eq. ( |D13| ) and Eq. ( |D14| )], and the form factor 

(Q) = 


(Dll) 


(D12) 


The final result for the structure factor, Eq. ( jDlOj ), is strik¬ 
ingly simple. We learn that, evaluated near to a reciprocal 
lattice vector Q, the structure factor S^{q) measures correla¬ 
tions of both P+ (q) and P_ (q), where each is projected onto 

^ X ^ X 

the vector Aq ^ and Aq _, respectively. These correlations 
are mixed with a weight controlled by the form factors F^{Q) 
and (Q). The form factors are, in turn, fully determined by 
the vectors Aq ^ and Aq _. These vectors therefore control 
both the “selection rules” which determine which of the fields 
P+ and P_ is manifest in the scattering around a given recip¬ 
rocal lattice vector, and the way in which these fields couple 
to a given experimental probe. 

All that remains is to substitute the appropria te Aq in the 
expression for S^{q) [Eq. (DIO)], using Eq. (Dll). In the 
specific case of the Ising structure factor 5'ising(q) [Eq. (|29|)], 
we have 

Af®(Q,q) 

3 3 

= exp(-iQ ■ r^) E I—exp[-i(Q + q) • df/2] 

i=0 


2 = 0 


± exp(-zQ • rc) E exp[-z(Q + q) • df /2] , 

(D13) 


Appendix E: Calculation of the dispersion of emergent photons 
on the ice Ih lattice 


The Hamiltonian of the U{1) gauge theory on the py- 


rochlore lattice is 


qj f 2 , 

^U(l) - ^ 2^ Kv' + ^ 

E 

(rr'>6CS 

(rr')eMS 

fC .-,2 F' 

F— ^ [Vo X A] + ^ 

E [Vo X (El) 

Oei 

Oeii 


where the sum 50(rr')ecs centre symmetric 

oxygen-oxygen bonds, 50(rr')eMS mirror symmet¬ 

ric oxygen-oxygen bonds [cf. Eig. [31, 
quettes normal to the optical axis and runs over pla- 

quettes parallel to the optical axis [cf. Eq. 0], We may 
condense this as 




(E2) 


where is a sum over primitive unit cells and the sums over 
m and p are over bond and plaquette midpoints in a single unit 
cell respectively. Eor the eight component objects a and /3 we 
have 


a = 

/9 = 


W W 
1,1,1,^,1,1,1,^ 


/r /c r r /c r \ 
V/c’/c’/c’^’/c’/c’/c’7' 


We will use the notation 


(E3) 

(E4) 


where the vectors d^, dp are given in Eqs. (A9)-(A15). 


Meanwhile, for the proton structure factor 


-E(r,m) — (E5) 

^(r,m) ~ ^r,r+dm (E6) 
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where is a bond vector in the direction of one of the eight 
inequivalent bonds in the unit cell. 


fields in terms of photon operators. This is 

A 








k A=1 




WA(k) 

(^exp [-ik ■ (r + d„/2)] r?„A(k)a;^(k) 

+ exp [zk • (r + d„,/2)] ?7^„(k)4(k)) (E7) 

T 


(r,m) 


= 


EE 

k A=1 


/^A(k) 

OLmJA 


^exp [—ik • (r + dm/2)] ? 7 mA(k)a;^(k) 

- exp [ik ■ (r + dm/2)] r?^m(k)aA(k)) (E8) 

fulfil the com- 


It is easy to check that that Eqs 
mutation relationship 


[-^(r,m): ^(r,m)] (E9) 

As in Ref. |[77]| , we may write Vq x A^m in terms of our 
Fourier decomposition as a matrix Z(k) acting on the vectors 

We need now to write down an decomposition of these 

I 

We have 


(Vo X ^)(r,p) — 



exp[-ik ■ (rp)]a^(k) ^/^Zpm{k)^/(^r]m\{k) 

m 


+ exp [zk • (rp)] ^(k) y] z?JmV^-^pm(E)*\/^ f (ElO) 


The matrix ^(k) is 

I(k) = 


0 

2isin(k-coi) 2i 

sin(k • C02) 

2i sin(k • cqs) 

0 

0 

0 

0 \ 

2i sin(k • Cio) 

0 

g-ik-Ci2 

_g—*k-Ci 3 

0 

0 

g—ik-Ci6 

_g-ik-Ci7 

2i sin(k • C20) 

_g-*k-C2i 

0 

g-ik-C 23 

0 

_g-*k-C 25 

0 

g-ik-C 27 

2i sin(k • C30) 

g-ik-Cai 

_g—ik-C 32 

0 

0 

g—ik-C35 

_g —*k-C 36 

0 

0 

0 

0 

0 

0 

2isin(k-C45) 2 isin(k-C 46 ) 2isin(k-C47) 

0 

0 

_g—ik-C 52 

g—ik-C53 

2i sin(k • C54) 

0 

_g ^k'C 56 

g—ik-C57 

0 

g-ik-Cei 

0 

_g—'ik-C 63 

2i sin(k • 054) 

g ''^k‘C 65 

0 

_g—ik-C 67 

V 0 

_g-ik-C7i 

g-ik-C 72 
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Where the vectors c^, Onm are vectors joining the central 
points of plaquettes to the bond midpoints around the outside. 
These vectors are defined in Appendix]^ along with the num¬ 
bering convention for sites and plaquettes. 

Defining 


(Ell) 


we may write 


v^(VoxA)„,„=7-5:5: 


u 


k A^ V 

x<^ exp[—zk- (rp)]a;^(k) x ^Ypm{k)r]m\{k) 

^ m 

+ exp [zk • (rp)] 4(k) x ?7)[m^Ap(k)|- (E12) 

m. ' 


— ‘sj'sjf^p 
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FIG. 11. (Color online). Identification of the excitations of the lattice gauge theory, ?^ u(i) [ Eg. ( |E1| )], with quantised fluctuations of the 
fields P+ and P-. (a) Dynamical structure factor of the uniform polarisation, P+ [Eq. ( |E25| )], calculated within the lattice gauge theory, 
(b) Equivalent results for the dynamical structure factor of the staggered polarisation, P_. Eor g ^ 0, the correlations of P+ are directly 
associated with the gapless, birefringent photons of the lattice gauge theory, while the correlations of P_ are associated with its gapped optical 
modes. This identification breaks down at shorter wavelengths, where the correlations of P+ and P- have contributions from all modes of 
the lattice gauge theory. Calculations were carried out for the symmetric choice of parameters, U = U' and JC = JC'. 


The matrix ^(k) is not guaranteed to be Hermitian, so we 
cannot necessarily form an orthonormal set of vectors r]^ (k) 
from its eigenvectors. However, the matrix ^^(k) • ^(k) is 

Hermitian and has positive semi-definite eigenvalues, so we 
may write 

V^(k) ■ V(k)^^(k) = Cx0^)\(k). (E13) 


Squaring and summing over r and p we have 


2 





k A 


^jj^|aA(k)aA(-k) + 4(k)4(-k) 

+aj^(k)aA(k)+aA(k)at^(k)|. (E14) 


Similarly 


^ am {Ermf = ^ E E 

(e,m) k A=1 


I - aA(k)aA(-k) - at^(k)4(-k) 
+aj^(k)aA(k) + aA(k)aJ^(k)|. (E15) 


Inserting this into the Hamiltonian [Eq. ( |E1[ )] results in 

•Hih = EE^^(^) f“A(k)aA(k) + (E16) 

where the dispersion co’A(k) is fixed by requiring the off- 
diagonal terms to vanish 

a;A(k) = v^|CA(k)|. (E17) 

The functions (x (k) are found by numerical diagonalization 
of Z^(k) • E(k). Eour of the modes A are unphysical, zero- 

energy modes (cf. the two unphysical zero energy modes 
which occur in the pyrochlore case (371). The remaining four 
modes are now non-degenerate, which is again in contrast to 
the cubic symmetry case. Two of these modes are gapless 
and linearly-dispersing, and therefore recognisable as pho¬ 
tons, while the other two modes are gapped and are associ¬ 
ated with quantised fluctuations of the classical field P_ (see 
Section|I^. The identification of the gapless modes with fluc¬ 
tuations of P+ and the gapped modes with fiuctuations of P_ 
is illustrated in Eig. im 

The time evolution of the electromagnetic fields is given 
by the time evolution of the photon operators aA(q),a^(q), 
which, since the photons are eigenstates of Hu(i )5 is simply 

4(q)(t) = aAq)(t) = 

(E18) 
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Therefore 


rocal lattice vector Q 


(E19) 

and the dynamical structure factor for the electric fields is 
= ^ H “ ^A(q))(l + nB{uj)) + 

A 

J(u;+CUA(q))nB(-^)] ^ ^ .^ VmXVXn • 

dmCiJA 

(E20) 


Since the electric field is in one-to-one correspondence 
with the proton configuration, the dynamical structure factor 
cc) also determines the scattering of neutrons or X- 
rays from protons in water ice. The relevant structure factor 
for the coherent scattering of neutrons is given by 1^ 


5'coh(q,^) 



(exp(—iq • Rj(0)) exp(iq • Rj/(t)) .) (E21) 


Within our treatment, the position of the proton on bond j of 
the lattice, Rj (f), is given by 


P±(q,^) = \^ exp{-iujt) X 
V J—oo 

exp(zQ • ta) X exp(-zq • r)P± (r^, t)+ 

reA 

exp(zQ • rc) X exp(-iq • r)P±(rc-, t) . 

rec 

(E25) 

where we have chosen a definition such that P±(q, oo) is the 
same in all Brillouin zones. The vectors r^, vc are the po¬ 
sitions of the A and C oxygen vertices within a primitive 
unit cell. Restricting the sum to the A and C sublattices of 
oxygen ions means that each bond only contributes to the 
sum once. In Eig.[^we plot (P+(—q, —co) • P+(q, cc)) and 

{P-(-q, -w) •P4q,w)). 

The results of this analysis are shown in Eig. At long 
wavelengths, the correlations of P+ are directly associated 
with the gapless emergent photons, while the fiuctuations of 
P_ are associated with the gapped optical modes. This con¬ 
firms the conclusions of the symmetry-based analysis given 
in Section JVB] 


Appendix F: Comparison with numerics 
1. Classical Monte Carlo simulation 


R-iW =r? + aj(Tj(^) 


(E22) 


where is the bond midpoint, describes the displacement 
of the two proton sites on the bond from the bond midpoint 
and (Tj{t) = zbl is the Ising variable describing the bond po¬ 
larisations. 


Using Eqs. ( |E21| )-( [E^ ) and the relationship between the 
Ising variables aj (t) and the electric fields Ej (t) 


(7j{t) = 2Ej{t) (E23) 

we find that the diffuse contribution to 5'coh (q? is given by 

5'TOh (q> ‘^) = 4 X ®“(q ■ ®“(q ■ "(q> A 

m,n 

(E24) 


where m, n index the eight sublattices of bonds and the dis¬ 
placement of a proton from the bond midpoint of bond m is 
±a^. 

We have also calculated the structure factor of the 
fields P+(r) and P-(r) within the lattice gauge theory 
^u(i) [Efi- (ET|]- To do this, by analogy with Eq. (|C14), 
we introduceme Eourier transform of the fields P±(qj for 
wavevector q = Q + q, in a Brillouin zone centred on recip- 


In this Appendix we compare our calculations of the corre¬ 
lation functions in the classical ice problem with the results of 
Monte Carlo simulations for a cluster of 1024 oxygen-oxygen 
bonds. 


To obtain numerical results for the correlations in the ab¬ 
sence of quantum tunnelling it is necessary to take an equally 
weighted average over ice rule configurations. This is accom¬ 
plished numerically by starting from a configuration with zero 
total electric polarisation and then acting randomly with the 
hexagonal ring exchange operators [cf. Eq. Q] to gener¬ 
ate new configurations. Since these operations preserve total 
electric polarisation the average obtained in this way only in¬ 
cludes states of vanishing total polarisation. However, this set 
of states is representative of the manifold as a whole, as in the 
case of spin ices. 


In Eig. [T^we compare results for the Ising structure factor 
5'ising(q) [Eq. (|29])] between these Monte Carlo calculations 
[(a)-(c)] and the projection method outlined in Appendix 
[(d)-(f)]. Noticeable differences are visible for q located ex¬ 
actly at Brillouin zone centers. This is caused by the restric¬ 
tion to states of vanishing total polarisation in the simulations. 
However, for any q not exactly at a Brillouin zone center 
there is very good agreement between the theory calculation 
and Monte Carlo, strongly validating our understanding of the 
classical ice problem. 
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FIG. 12. (Color online). Ising structure factor Arising(q) [Eq. j^ ] for a classical model of ice Ih, calculated from classical Monte Carlo 
simulation [(a)-(c)] and from lattice theory [(d)-(f)]. The theory calculation is performed using the method outlined in Appendix [C] which is 
based on the method described for spin ice in Ref. l54l . This method consists in writing the ice rule constraints as orthogonality conditions 
in Fourier space. In the long wavelength limit these conditions reduce to those obtained from the continuum field-theory presented in Section 
|III A| [Eqs. (p^-|T8|]. Near the zone centres the correlations are well decribed by a combination of pinch point singularities, reflecting the 
algebraic correlations of P+(q) and smooth features reflecting the short ranged correlations of P- (q). 


2. Quantum Monte Carlo simulation 


In order to validate our description of quantum ice Ih, we 
have compared the predictions of the lattice gauge theory 
^u(i) [Eq. (66)] for equal-time correlation functions, with 
the results of variational quantum Monte Carlo (VMC) simu- 


hexagonal 

tunnelling 


[Eq. 4 ] 


lations of the microscopic model 

VMC simulations were carried out for an ofEfiorhombic 
cluster of 1024 oxygen-oxygen bonds with a one-parameter 
variational wavefunction of the form 


1?/)^^'') = exp[aA//]|VAo) (FI) 

where the operator A// measures the number “flippahle” 
plaquettes (of both type / and II), a is a variational pa¬ 
rameters, and I'lpo) is an equal weight superposition of all 
ice configurations within a given flux-sector of the Hilbert 
space. This variational wave funtion correctly describes the 
physics of quantum liquids based on a lattice U(l) gauge the¬ 
ory l[37l[39ll43ll44l . but is heavily biased toward liquid, rather 


than ordered ground states ll44l . Since our intention here is 
to investigate the properties of a proton-liquid in water ice, 
rather than to investigate possible ordered ground states, the 
wave function Eq. is sufficient for our purposes. 

In Fig. [^we show the equal-time “Ising” structure factor 

S,,{q,t = 0) = (S^(-q,t = 0)S^(q,t = 0)) (F2) 


in the (/i, /c, 0) plane of reciprocal space. Results were calcu¬ 
lated from the lattice gauge theory 'Hu(i) [Eq. (66)], for pa¬ 
rameters 


K = K' = IC = IC' = 1 (F3) 

and from the microscopic model [Eq. (|^], for pa¬ 

rameters 

gi= 92 = 1 (F4) 

using VMC, as described in (221111. A fleld-renormalization 
of S^, allowed in principle within the lattice gauge theory |[36| 
ITtI . has also been set equal to one. 
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FIG. 13. (Color online). Comparison between the predictions of the lattice gauge theory ?^u(i) [Eq. ([66|], and variational Monte Carlo 
simulation of [Eq. 1^] for the “Ising” structure-factor S'ising(q, t = 0) [Eq. jp^], for the parameter sets given in Eq. ( [f 3| and 

Eq. There is excellent agreement between the two methods, validating our description of the proton-liquid in water ice Ih in terms of a 
quantum U(l) lattice gauge theory. 


The agreement between the results of the two methods is 
excellent, confirming that the lattice gauge theory correctly 
describes the liquid phase of the microscopic model. 


The elastic contribution is simply 

S'inc.ei (q, w) = (5(a;) cos^ (q • aj) (G6) 


Appendix G: Calculation of the incoherent scattering cross 
section 

The structure factor for inelastic, incoherent, neutron scat¬ 
tering from a set of nuclei at located at positions Rj is ll99l 

1 

5inc(q,a;) = —dte"-* 

(exp(—iq • Rj(0)) exp(iq • Rj(f))) (Gl) 

We approximate the position of the proton on bond j of the 
lattice, Rj(f), to be given by 

Rj(t) = (G2) 

where is the bond midpoint, aj describes the displacement 
of the two proton sites on the bond from the bond midpoint 
and = ±1 is the Ising variable describing the bond po¬ 
larisations. 

Since dj (t) = ±1 and r? and aj are constants we can write 


while the inelastic contribution is 

5'inc,inei(q,w) = y^sin^(q • aj){aj{-u)aj{Lo)) (G7) 

j 


For the purposes of comparison with experiments we need 
to integrate the momentum transfer dependence over angle, 
which gives 


j ^ 


^ sin(2Q|aj|) \ 

2Q|aj| / 


sin(2Q|aj|) ^ 
2Q\a,\ ) ^ 


(G8) 


-w)crj(w)) 


(G9) 


where Q = |q|. 

The local correlation function {aj{—uj)(Tj{uj)) can be 
rewritten as a sum over Fourier space, so we have 

8 


exp(iq • Rj (t)) 

= exp(iq • r^) [cos(q • a^) + 'io-j{t) sin(q • a^)] 


Since 


{aj{t)) = 0 


corresponding to the elastic and inelastic contributions to the 
incoherent scattering 

^inc(q5^) ~ ^inc,el(q5^) ^inc,inel (q? (GS) 


i=l 


S'inc.inei(q, sin^(q • a*) y^(cri(-q', , u)) 

q' 

(GIO) 




(G3) 

QPOW 

^inCjinel 


(Q,w) = 

(G4) 




(Gll) 

where the sum over i now runs over the eight sublattices of 
bonds. 
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Within the lattice gauge theory description the correlations 
of the Ising variables a is directly related to the correlations 
of the electric field E 

(o-i(-q',-w)(Tj(q',w)) = 4(£'i(-q',-w)£'i(q',w)) 

(G12) 

SO calculating the correlations of E in the gauge theory en¬ 
ables us to calculate the incoherent scattering from protons, 
as shown in Fig. The effects of finite temperature, which 
lead to the thermal excitation of photons, can also be included, 
as described in Oil. 

Calculations were carried out for the symmetric parameter 
set given in 


U = U' , }C = JC^ 

[cf. Eq. ([68])], with energy scale 

Vwc = 0.018 meV , 

[cf. Eq. ( [8^ ] . The one remaining free parameter of the theory, 
is used to fix the average normalisation of the electric 


fields of the gauge theory [cf. Eq. (|65ll]. At T = 5K this gives 

= 47.8. (G13) 

Experimental measurements are carried out with finite en¬ 
ergy resolution. To make a comparison with experiment, it 
is therefore necessary to convolute both the elastic response 
(a delta-function in energy) and the prediction for incoherent 
inelastic scattering with a representation of the experimental 
line-shape. We adopt the simplest representative line shape, a 
Gaussian 



5Xp(^) 




(G14) 


with energy width 


S = 0.03 meV , 

chosen as representative of the elastic line shown in Ref. [T8| 
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